Robust Control of a Multivariable Experimental Four-Tank System
The goal is to explore some results of the two papers
- K. H. Johansson. The Quadruple-Tank Process.- A multivariable laboratory process with an adjustable zero. IEEE Control Systems Technology, 2000;
- Rajanikanth Vadigepalli, Edward P. Gatzke, and Francis J. Doyle III, Robust Control of a Multivariable Experimental Four-Tank System. Ind. Eng. Chem. Res. 2001.
regarding the analysis of the properties of a multivariable Four-Tank System with 2 inputs and 2 outputs. We will also provide controllers for robust stabilizing the linear tangent model of the non linear system around the origin. Our dissertation will follow these keypoints:
- analysis of the eventual Minimum Phase property of the linearized system depending on the choice of parameters;
- finding the directions of the zeros and the poles of the linearized system and explore possible implications;
- showing the limitations due to the presence of eventual Non Minimum Phase zeros;
- introducing increasing uncertainty and trying to develop some robust decentralized controllers;
- exploring some decoupling approaches.
The goal is to control the levels
and
in the two lower tanks using the two pumps. The plants control inputs are the two input voltages of the pumps,
and
. The fòow produced by pump i corresponding to a voltage
is
. The i-th valve distributes the output ow of the pump by sending
to the lower tank and the remaining
to the upper tank. Clearly
. The gravity acceleration is g while
is the cross section area of tank i and
is the cross-section area of the outlet hole of each tank. A sensor provides an output voltage proportional (with constant
) to the water height
, therefore the two measured outputs (two sensors with same
) are
Denoting with
i the height at the equilibrium in each tank, the time costants are given by: 
The values of these parameters are taken from the Johansson' paper.
1. Poles and Zeros location and Minimum Phase property analysis.
The control object of our analysis is the MIMO transfer futnction matrix of the linearized model around an equilibrium point:
Let us first of all compute the poles and the zeros.
We know that in a MIMO system the tranfer function's poles can be computed as the roots of the pole polynomial. The pole polynomial
, i.e. the least common denominator of all non-zero minors of all orders of G(s). The minors of order 1 are the entries
of the matrix, while the only minor of order 2 is the transfer function determinant:
.so 
. Then we can state that the poles of the transfer function are
each of them with multeplicity equal to 1.
Similar considerations can be done for the computation of the zero. In the general case, the zeros can be computed as the roots of the zero polynomial
i.e. the greatest common divisor of all the numerators of all order r (normal rank) minors of G(s), provided that these have been adjusted so to have the pole polynomial as denominator. So if the normal rank of
is 2, we need to check only the minor of order 2, i. e. the determinant of the matrix. Additionally, becuse of the fact that the transfer function is square (p = q = 2), we could directly say that the solutions of
0 gives the transmission zeros which do not coincide with a pole. Let's check when the determinant is different from zero. If it is equal to zero only at few values (normal rank = 2), we can state that such values are the transmission zeros of the transfer function.
so:
a polynomial of order 2 whoose roots are:
and consequently the term
=
(
-
) So, for the physics of the system, parameters will never be such that these two solutions are complex, and they costitutes the only two values with respect to which the matrix G(s) loses rank. So G(s) has normal rank equal to 2, and that solutions are the symbolical two trasmission zeros of the trasfer function.
Let us validate our computations using the calculator:
syms g1 g2 k1 k2 n1 n2 t1 t2 t3 t4 s A1 A2 A3 A4 kc % inizializing the symbolic parameters of the system
G = [(g1*k1*n1/(1+t1*s)), ((1-g2)*k2*n1)/ ((1+t3*s)*(1+t1*s)); ((1-g1)*k1*n2)/ ((1+t2*s)*(1+t4*s)), (g2*k2*n2)/(1+t2*s)] % inizializing the plant
G =

d = det(G)
d =

zi = solve(eqn, s, 'ReturnConditions', true);
simplify(zi.s)
ans =

These are the result we obtained, barring semplifications.
From the above analysis, it is easy to verify the properties of those zeros.
- First of all, from the expresion of
, we can see that if
c is always equal to 0, so
and
is Re[ ] < 0. - secondly, if
the quantity -4ac is always >0, because the quantity
is greater that
, so there is a positive zero and a negative one (system is non minimum-phase). - thirdly, if
the quantity -4ac is always < 0, because the quantity
is smaller that
, so there are two negative zeros (system is minimum-phase). Naturally, since
we can have at most that
So the minimum-phase conditions is simply
< 2.
Thinking about the physical interpetation these results make sense. We have a non minimum phase system if
, i.e. if the entering flow into the two lower thank ( the ones with respect to which we want to regulate the water level) is smaller than the total one entering the two upper tanks. Viceversa, we have minimum phase property is the upper tanks receive more flow tha the others. If instead the total flow going to the left tanks 1 and 3 is equal to the total flow going to the right tanks 2 and 4, we have to manage with some more difficulties derived by the presence of the zero in 0.
Note that:
.So, in the case in which
and
, we have that
and 
, so the two zeros of the transfer function are close to be coincident with the two poles
,
.
Note also that if one of the two or both
we have not transimition zeros in the transfer function.
2. Directions of the zeros and the poles, blocking property of the zeros.
In order to compute the zeros directions, we need a minimal realization {A,B,C,D}: without it, we could not use the useful results about the system matrix.
Let us check if the proposed realization in the Johansson's paper is minimal.
To do this, we need to check the observability and the controllability properties of the natural modes of that system. If the system is fully controllable and fully observable, we can state that the realizazion {A,B,C,D} is minimal.
So let us compute the Observability and the Controllability Matrices and check if they are full rank.
A = [ -1/t1, 0 , A3/(A1*t3), 0; 0 , -1/t2, 0 A4/(A2*t4); 0, 0, -1/t3, 0; 0,0,0, -1/t4]; % matrix A
C = [kc, 0,0,0; 0, kc,0,0]; % matrix C
B = [(g1*k1)/A1, 0; 0, (g2*k2)/A2; 0, (1-g2)*k2/A3; (1-g1)*k1/A4, 0]; % matrix B
Obv = [C; C*A; C*A_2; C*A_3]; % observability matrix
Con = [B , A*B, A_2*B, A_3*B]; % controllability matrix
Both matrices have rank = 4, equal to the dimention of the state. This implies that the system is fully controllable and observable. The state space representation is then a minimal realization of the transfer function G(s).
Let us now substitute the symbolical quantities with the values used in the Johansson paper and try to compute these quantities.
Our considerations on the zeros and poles directions will regard the operating points
and
and around which the system exhibits respectively minimum-phase (MP) and non-minimum-phase (NMP) characteristics.
t1_mp = (A1_n/c1_n) * sqrt(2*h01_mp/g_n); % s
t2_mp = (A2_n/c2_n) * sqrt(2*h02_mp/g_n); % s
t3_mp = (A3_n/c3_n) * sqrt(2*h03_mp/g_n); % s
t4_mp = (A4_n/c4_n) * sqrt(2*h04_mp/g_n); % s
n1_mp = (kc_n* t1_mp)/A1_n;
n2_mp = (kc_n* t2_mp)/A2_n;
t1_nmp = (A1_n/c1_n) * sqrt(2*h01_nmp/g_n); % s
t2_nmp = (A2_n/c2_n) * sqrt(2*h02_nmp/g_n); % s
t3_nmp = (A3_n/c3_n) * sqrt(2*h03_nmp/g_n); % s
t4_nmp = (A4_n/c4_n) * sqrt(2*h04_nmp/g_n); % s
n1_nmp = (kc_n* t1_nmp)/A1_n;
n2_nmp = (kc_n* t2_nmp)/A2_n;
% MP state space inizialization
A_mp = eval(subs(A, [A1, A2, A3,A4, t1, t2, t3 ,t4], [A1_n, A2_n, A3_n, A4_n, t1_mp, t2_mp, t3_mp, t4_mp]));
B_mp = eval(subs(B, [A1, A2, A3,A4, g1, g2, k1 , k2], [A1_n, A2_n, A3_n, A4_n, g1_mp, g2_mp, k1_mp, k2_mp]));
C_mp = eval(subs(C, kc, kc_n));
G_sys_mp = ss(A_mp, B_mp, C_mp, D_mp);
% NMP state space inizialization
A_nmp = eval(subs(A, [A1, A2, A3,A4, t1, t2, t3 ,t4], [A1_n, A2_n, A3_n, A4_n, t1_nmp, t2_nmp, t3_nmp, t4_nmp]));
B_nmp = eval(subs(B, [A1, A2, A3,A4, g1, g2, k1 , k2], [A1_n, A2_n, A3_n, A4_n, g1_nmp, g2_nmp, k1_nmp, k2_nmp]));
C_nmp = eval(subs(C, kc, kc_n));
G_sys_nmp = ss(A_nmp, B_nmp, C_nmp, D_nmp);
% inizializing the transfer function of the plant with MP property
Gs_mp = [(g1_mp*k1_mp*n1_mp/(1+t1_mp*s)), ((1-g2_mp)*k2_mp*n1_mp)/ ((1+t3_mp*s)*(1+t1_mp*s)); ((1-g1_mp)*k1_mp*n2_mp)/ ((1+t2_mp*s)*(1+t4_mp*s)), (g2_mp*k2_mp*n2_mp)/(1+t2_mp*s)]
Gs_mp =
From input 1 to output...
2.61
1: ----------
62.7 s + 1
1.41
2: ----------------------
2709 s^2 + 120.3 s + 1
From input 2 to output...
1.5
1: ----------------------
1498 s^2 + 86.59 s + 1
2.837
2: -----------
90.34 s + 1
Continuous-time transfer function.
Model Properties
% inizializing the transfer function of the plant with NMP property
Gs_nmp = [(g1_nmp*k1_nmp*n1_nmp/(1+t1_nmp*s)), ((1-g2_nmp)*k2_nmp*n1_nmp)/ ((1+t3_nmp*s)*(1+t1_nmp*s)); ((1-g1_nmp)*k1_nmp*n2_nmp)/ ((1+t2_nmp*s)*(1+t4_nmp*s)), (g2_nmp*k2_nmp*n2_nmp)/(1+t2_nmp*s)]
Gs_nmp =
From input 1 to output...
1.524
1: -----------
63.21 s + 1
2.556
2: ----------------------
5128 s^2 + 147.5 s + 1
From input 2 to output...
2.451
1: ----------------------
2466 s^2 + 102.2 s + 1
1.597
2: ----------
91.4 s + 1
Continuous-time transfer function.
Model Properties
Let's first of all compute the zeros of the system in the way explained in point 1, both for MP and NMP cases, i.d. seeing the values for wihich 
det_mp = (k1_mp*k2_mp*n1_mp*n2_mp*(g1_mp + g2_mp + g1_mp*g2_mp*s*t3_mp + g1_mp*g2_mp*s*t4_mp + g1_mp*g2_mp*s^2*t3_mp*t4_mp - 1))/((s*t1_mp + 1)*(s*t2_mp + 1)*(s*t3_mp + 1)*(s*t4_mp + 1));
zi_mp = zero(det_mp) % computation of the 2 minimum phase transmission zeros
det_nmp = (k1_nmp*k2_nmp*n1_nmp*n2_nmp*(g1_nmp + g2_nmp + g1_nmp*g2_nmp*s*t3_nmp + g1_nmp*g2_nmp*s*t4_nmp + g1_nmp*g2_nmp*s^2*t3_nmp*t4_nmp - 1))/((s*t1_nmp + 1)*(s*t2_nmp + 1)*(s*t3_nmp + 1)*(s*t4_nmp + 1));
zi_nmp = zero(det_nmp) % computation of the 1 minimum phase and the 1 non minimum phase transmission zeros
As we can see, the signs of the zeros are coherent with what we expected. Let's now compute the system matrix
and see if its zeros coincides with the zero of the transfer function. We have to remember that since the system is square (p = q), left and right zeros coincides. Furthermore, since the realization is minimal, the zeros of the system matrix
with respect to which it drops rank are ONLY the transmition zeros, there is no presence of invariant zeros. sys_matrix = [s*eye(4)- A, -B; C, D] % inzializing the system matrix
sys_matrix =

sys_matrix_mp = eval([s*eye(4)-A_mp, -B_mp; C_mp, D_mp]); % inzializing the system matrix MP
sys_matrix_nmp = eval([s*eye(4)-A_nmp, -B_nmp; C_nmp, D_nmp]); % inzializing the system matrix NMP
det_sys_mat_mp = eval(det(sys_matrix_mp)); % determinant of MP system matrix
det_sys_mat_nmp = eval(det(sys_matrix_nmp)); % determinant of NMP system matrix
eq5 = det_sys_mat_mp == 0; % finding values with respect to which the system matrix drops rank
eq6 = det_sys_mat_nmp == 0;
zi_mp= eval(solve(eq5,s)) % 2 MP zeros
zi_nmp= eval(solve(eq6,s)) % 1 MP 1 NMP zeros
Which proves our expectations. Let us further validate our results using the MATLAB function tzero() to see if there is matching of computations.
Verified that the computed two pair of zeros are right, we can now find the zero directions solving the two problems for every zero:
where
is the system matrix computed at one possible zero,
,
, respectively the right and left generalized eigenvector associated to z. The vector
costitutes the zero state direction ,
the zero input direction,
the zero output direction. We can reformulate this two problems as a generalized eigenvalue problem:
with
and
, and find
easily. Let us follow this approach.
M = [eye(4), B_mp*0;C_mp*0, D_mp]; % inizialization of M
G_mp = [A_mp,B_mp; C_mp, D_mp]; % inizialization of G (MP zeros)
G_nmp = [A_nmp,B_nmp;C_nmp, D_nmp]; % inizialization of G (1 NMP zeros)
[na,za,ma] = eig(G_mp,M) % computing n and m
0 0 0 0 0.0000 0
0 0 0 0 0 0.0000
0 1.0000 -1.0000 -0.0000 -0.0000 0.0000
0 -0.9714 -0.6361 -0.0000 0.0000 0
1.0000 -0.5028 0.5028 0 1.0000 -0.0395
0 0.5156 0.3377 -1.0000 -0.0000 1.0000
Inf 0 0 0 0 0
0 -0.0172 0 0 0 0
0 0 -0.0580 0 0 0
0 0 0 Inf 0 0
0 0 0 0 Inf 0
0 0 0 0 0 Inf
0.0000 0.3750 0.3750 0 0 0
0.0000 -0.4847 0.7401 -0.0000 0 0
0 0.6361 -0.9714 0 0 0
0.0000 -1.0000 -1.0000 0 0 0
1.0000 -0.0009 -0.0316 0 1.0000 0
0.0409 0.0059 -0.0695 -1.0000 0 1.0000
Here we can see the generalized eigenvectors associated to the two zeros
and
for the equilibrium point
. The columns 2,3 of the matrices
,
are respectively the two right and left generalized eigenvectos; while the first 4 elements of the 2° and 3° columns of
are the zeros state direction
and the last two elements of the 2° and 3° columns of
and
are respectively the zero input and output directions
,
. The following are the same quantities for the equilibrium point 
[nb,zb,mb] = eig(G_nmp,M)
0 0 0 0 0.0000 0
0 0 0 0 0 0.0000
0 1.0000 -1.0000 0.0000 0 -0.0000
0 -0.9716 -0.7740 0 0.0000 0.0000
1.0000 -0.5316 0.5316 -0.0557 1.0000 0.0946
0 0.4953 0.3946 -1.0000 -0.0000 1.0000
Inf 0 0 0 0 0
0 0.0128 0 0 0 0
0 0 -0.0562 0 0 0
0 0 0 Inf 0 0
0 0 0 0 Inf 0
0 0 0 0 0 Inf
0.0000 0.6755 0.5381 0 0 0
0.0000 -1.0000 1.0000 -0.0000 0 0
0.0000 0.4508 -0.4508 0 0 0
0 -0.5824 -0.4639 0 0 0
1.0000 0.0386 -0.0435 0 1.0000 0
0.0481 -0.0474 -0.0906 -1.0000 0 1.0000
Let us now extract the velues specified before.
csi_z_mp = na(1:4, 2:3) % zero state directions for the two MP zeros of P-
0 0
0 0
1.0000 -1.0000
-0.9714 -0.6361
v_z_mp = na(5:6, 2:3) % zero input directions for the two MP zeros of P-
-0.5028 0.5028
0.5156 0.3377
csi_z_nmp = nb(1:4, 2:3) % zero input state directions for the two zeros (one MP one NMP) of P+
0 0
0 0
1.0000 -1.0000
-0.9716 -0.7740
v_z_nmp = nb(5:6, 2:3) % zero input directions for the two zeros (one MP one NMP) of P+
-0.5316 0.5316
0.4953 0.3946
w_z_mp = ma(5:6, 2:3) % zero output directions for the two MP zeros of P-
-0.0009 -0.0316
0.0059 -0.0695
w_z_nmp = mb(5:6, 2:3) % zero output directions for the two zeros (one MP one NMP) of P+
0.0386 -0.0435
-0.0474 -0.0906
eta_z_nmp = mb(1:4, 2:3) % zero output state directions for the two zeros (one MP one NMP) of P+
0.6755 0.5381
-1.0000 1.0000
0.4508 -0.4508
-0.5824 -0.4639
From the inspection of the component of
, we can see that the direction of the NMP zeros is not a unit vector, so the effect of the RHP zero is distributed between both outputs.
Now we want to show the blocking property of these two pairs of transmition zeros.
In other words, we want to show that the outputs of the two system
exept for some small errors due to numerical digit approximations (always present in Matlab), provided that the inputs are set as
and
and the inizial states
. This means that every couple of inputs with such directions
, with the form of an exponential with grow rate equal to
, will not be reproduced in non of the outputs provided that the system starts to evolve from the zero state direction
, thanks to the presence of the associated transmition zeros
.
From the inspection of the component of
, we can state both for the MP and the NMP case that, in order to the exhibition by the system of the blocking property of the zeros, the distances of the levels of the first two tanks
and
from the height at the two equilibrium point
and
have to be equal to zero, the remaining 3 and 4 have to differ by the quantity specified in the other 2 elements of
, and the difference between the current voltages provided and the ones at the two equilibria has to be a signal with exponential form, grow rate equal to
, and the proportions in which the exponential should be present at the corresponding input equal to the components of
. Let us try to visalize all of these things by simulations.
t = 1:0.001:1000; % inizializing the time span
u1_mp = [zeros(length(t),1), zeros(length(t),1)]; % inizializing u_1 at P_
u2_mp = [zeros(length(t),1), zeros(length(t),1)]; % inizializing u_2 at P_
u1_nmp = [zeros(length(t),1), zeros(length(t),1)]; % inizializing u_1 at P+
u2_nmp = [zeros(length(t),1), zeros(length(t),1)]; % inizializing u_2 at P+
% building the exponential signals with directions v_z :
u1_mp(i,1) = exp(zi_mp(2)*t(i))* v_z_mp(1,1);
u1_mp(i,2) = exp(zi_mp(2)*t(i))* v_z_mp(2,1);
u2_mp(i,1) = exp(zi_mp(1)*t(i))* v_z_mp(1,2);
u2_mp(i,2) = exp(zi_mp(1)*t(i))* v_z_mp(2,2);
u1_nmp(i,1) = exp(zi_nmp(2)*t(i))* v_z_nmp(1,1);
u1_nmp(i,2) = exp(zi_nmp(2)*t(i))* v_z_nmp(2,1);
u2_nmp(i,1) = exp(zi_nmp(1)*t(i))* v_z_nmp(1,2);
u2_nmp(i,2) = exp(zi_nmp(1)*t(i))* v_z_nmp(2,2);
% simulating from inizial conditions csi_z:
y1 = lsim(G_sys_mp, u1_mp,t, csi_z_mp(:,1));
y2 = lsim(G_sys_mp, u2_mp,t, csi_z_mp(:,2));
y3 = lsim(G_sys_nmp, u1_nmp,t, csi_z_nmp(:,1));
y4 = lsim(G_sys_nmp, u2_nmp,t, csi_z_nmp(:,2));
plot(t, y1, 'LineWidth',2); grid
title("MP case: blocking property of the zero z1 = -0.0172");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
plot(t, y2, 'LineWidth',2); grid
title("MP case: blocking property of the zero z2 = -0.0580");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','East')
plot(t, y3, 'LineWidth',2); grid
title("NMP case : blocking property of the zero z1 = -0.0562");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
plot(t, y4, 'LineWidth',2); grid
title("NMP case : blocking property of the zero z2 = 0.0128");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','NorthEast')
As we expected, the outputs are identically = 0 except for some small errors of the order of
due to the numerical approximations. In order investigate more in which way the effect of the RPH zero is distributed between the two outputs, we can also compute symbolically the zero output directions through the transfer function, as done in the Johansson's paper.
We start putting : 
res = wz'*G
res =

Then we solve for
as suggested in the Johansson's paper. eq7 = res(2) == 0
eq7 =

solution = solve(eq7,g2) % solving for gamma2
solution =

Then we substitute this values in the original expression, we solve for
we simplify for
. res1= eval(subs(res, g2, solution)) % substituting gamma2 in the original expression
res2 = collect(res1, wz1) % regrouping terms
s2 = solve(eq8, wz1) % solving for w_z1
s2 =

eq9 = wz1/wz2 == s2/wz2 % obtaining the Johansson's relation
eq9 =

This relation shows that as
decreases,
increases and the effect of the zero is more relevant in the first output. Viceversa, as
increases
decreases and the zero has a more relevant effect for
This is intuitive: as much we provide less flow to one tank, harder to control is its level (for the increased effect of the NMP zero).
Let us now find the symbolical pole directions, remembering that
are the input poles directions and
the output poles directions, where
,
are respectively the left and the right eigenvectors associated with the eigenvalues of the matrix A. [U, Di] = eig(A) % computation of right eigenvectors
U =

Di =

[V, Di] = eig(A.') % computation of left eigenvectors (computed as columns)
V =

Di =

input_pole_vec = (V.')*B % each row is the input direction associate to i-th eigenvalue (i-th pole)
input_pole_vec =

output_pole_vec = C*U % each column is the output direction associate to i-th eigenvalue (i-th pole)
output_pole_vec =

It is interesting to evaluate these quantities in the parameters combination we analyzed in the point 1, the one in which
We saw that in such a case we have two zero-pole coincidences in 
We know that such a coincidence will bring to a zero-pole cancellation conditionally to several condition.
In particular we have:
- a pole/zero cancellation with loss of observability iff
where
is the pole right eigenvector; - zero/pole cancellation with loss of controllability iff
where
is the pole left eigenvector.
So let us evaluate the system numerically with
and compute the poles and the zeros directions.. A_canc= eval(subs(A, [A1, A2, A3,A4, t1, t2, t3 ,t4], [A1_n, A2_n, A3_n, A4_n, t1_mp, t2_mp, t3_mp, t4_mp]));
B_canc = eval(subs(B, [A1, A2, A3,A4, g1, g2, k1 , k2], [A1_n, A2_n, A3_n, A4_n, 1, 1, k1_mp, k2_mp]));
C_canc = eval(subs(C, kc, kc_n));
[U_canc, D] = eig(A_canc) % computation of right eigenvectors
1.0000 0 -0.8503 0
0 1.0000 0 -0.8315
0 0 0.5263 0
0 0 0 0.5555
-0.0159 0 0 0
0 -0.0111 0 0
0 0 -0.0419 0
0 0 0 -0.0333
[V_canc, D] = eig(A_canc.') % computation of left eigenvectors
0 0.5263 0 0
0 0 0 0.5555
1.0000 0.8503 0 0
0 0 1.0000 0.8315
-0.0419 0 0 0
0 -0.0159 0 0
0 0 -0.0333 0
0 0 0 -0.0111
input_pole_vec_canc = (V_canc.')*B_canc % computation of input poles directions
0 0
0.0626 0
0 0
0 0.0581
output_pole_vec_canc = C_canc*U_canc % computation of output poles directions
0.5000 0 -0.4251 0
0 0.5000 0 -0.4158
Gs_canc = [(k1_mp*n1_mp/(1+t1_mp*s)), 0; 0, (k2_mp*n2_mp)/(1+t2_mp*s)] % transfer function
Gs_canc =

G_ss_canc = ss(A_canc, B_canc, C_canc, D_canc);
As we expected we have two zeros/pole coincidences in
. G_canc = [A_canc,B_canc; C_canc, D_canc]; % inizialization of G (gamma1 = gamma2 = 1)
[n_canc,z_canc,m_canc] = eig(G_canc,M) % computing n and m
0 0 0 0 0.0000 0
0 0 0 0 0 0.0000
0 0 1.0000 0 0 0
0 0 0 1.0000 0 0
1.0000 0 -0.3520 0 1.0000 0
0 1.0000 0 -0.3185 0 1.0000
Inf 0 0 0 0 0
0 Inf 0 0 0 0
0 0 -0.0419 0 0 0
0 0 0 -0.0333 0 0
0 0 0 0 Inf 0
0 0 0 0 0 Inf
0.0000 0 0 0 0 0
0 0.0000 0 0 0 0
0 0 1.0000 0 0 0
0 0 0 1.0000 0 0
1.0000 0 0 0 1.0000 0
0 1.0000 0 0 0 1.0000
v_z_canc = n_canc(5:6, 3:4) % zero input directions (gamma1 = gamma2 = 1)
csi_z_nmp = n_canc(1:4, 3:4) % zero state input directions (gamma1 = gamma2 = 1)
eta_z_canc = m_canc(1:4, 3:4) % zero state output directions (gamma1 = gamma2 = 1)
w_z_canc = m_canc(5:6, 3:4) % zero output directions (gamma1 = gamma2 = 1)
As we can see,
, so with dont's have a pole/zero cancellation with loss of observability. Additionally,
, so we have not a zero/pole cancellation with loss of controllability. So this proves that altought there are two zero/pole coincidences as , there are never cancellations and hidden dynamics in the process. This result is coherent with the fact that the realization {A,B,C,D} used is minimal.
3. Limitations due to the presence of NMP zeros.
The first limitation we want to explore is the fact that if G(s) has a NMP-zero at z with output direction
, then for internal stability of the closed-loop system it must be guaranteed the validity of the two following interpolation constraints : i. e. T must have a NMP-zero in the same direction as G(s), and S(z) has an eigenvalue of 1 corresponding to the left eigenvector
. Let's compute a simple stabilizing controller using
norm approach, withouth taking into account other things like specifications on the control effort or decoupling aims. A diagonal weight seems to be sufficient.
wp1 = (s/M1+ omb1)/ (s+omb1*A1);
wp = (s/M1+ omb1)/ (s+omb1*A1)* eye(2); % diagonal weight
[KHinf1, Ghinf1, gopt1] = mixsyn(Gs_nmp, wp, [], []);
Cloop1 = loopsens(Gs_nmp, KHinf1);
r_unit = [ones(length(t),1),ones(length(t),1)]; % dummy reference to check if the controller works
ycl1 = lsim(Cloop1.To, r_unit, t);
plot(t, ycl1, 'LineWidth', 2), grid
title("System's steady state with simple H_{inf} controller");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','NorthEast')
The closed loop system is stable. We can see that the interpolations conditions hold.
%T = (Gs_nmp)/(eye(2)+Gs_nmp); % inizializing complementary sensitivity matrix
S = Cloop1.So; % inizializing sensitivity matrix
T_z = evalfr(T, zi_nmp(2)); % computing at z = 0.128
S_z = evalfr(S, zi_nmp(2)); % computing at z = 0.128
w_z_nmp(:,1)' * T_z % showing the limitation
w_z_nmp(:,1)' * S_z % showing the limitation
We expect then a NMP zero in
in T(s) exhibiting as reasonable the blocking property, exactly like G(s). So let us solve the generalized eigevectors problem in order to compute the zeros input and output directions, and the zero state directions.
GT = [T_sys.A, T_sys.B; T_sys.C, T_sys.D]; % G matrix for T(s)
MT = [eye(14), T_sys.B*0; T_sys.C*0, T_sys.D*0]; % M matrix for T(s)
[nT,zT,mT] = eig(GT,MT)
0 0.0000 0.0000 0.0000 0.3665 -0.0798 0.1505 0.0000 0.0000 0.0000 -0.2490 -0.0000 -0.0623 -0.0000 -0.0000 -0.0000
0 0 -0.0000 -0.0000 0.9575 -0.0804 0.0000 -0.7002 0.0000 0.7002 0.0000 -0.0000 0.0000 -0.0009 0.0000 0.0000
0 0 0 -0.0000 -0.2661 -0.0983 -0.0000 1.0000 -0.0000 -1.0000 -0.0000 0.0000 -0.0000 0.0013 0.0000 0.0000
0 0 0.0000 0.0000 1.0000 0.0495 0.1155 -0.0000 -0.0000 0.0000 -0.1911 -0.0000 -0.0478 -0.0000 -0.0000 -0.0000
0 0 0 0.0000 -0.2779 0.0605 -0.1141 -0.0000 -0.0000 -0.0000 0.1887 0.0000 0.0472 0.0000 -0.0000 -0.0000
0 0 0.0000 0.0000 0.2428 0.0897 0.0000 -0.9125 0.0000 0.9125 0.0000 -0.0000 0.0000 -0.0012 -0.0000 -0.0000
0 -0.0000 -0.0000 -0.0000 0.0008 -0.7790 1.0000 0.0000 1.0000 0.3740 1.0000 0.6955 1.0000 -0.7131 0.0000 0.0000
0 -0.0000 -0.0000 0.0000 0.0003 1.0000 -0.9740 -0.0000 -0.8936 -0.5244 -0.9740 -1.0000 -0.9740 1.0000 -0.0000 -0.0000
0 -0.0000 -0.0000 -0.0000 -0.3665 0.0798 -0.0025 -0.0000 0.0000 -0.0000 -0.1820 0.0000 0.0549 -0.0000 0.0000 0.0000
0 0.0000 0.0000 0.0000 -0.9575 0.0804 -0.0000 0.7002 -0.0000 0.5802 -0.0000 0.1943 -0.0000 -0.0328 -0.0000 -0.0000
Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 -0.0562 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0.0128 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 -0.0158 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 -0.0109 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 -0.0256 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 -0.0109 0 0 0 0 0 0
0.0000 -0.0000 -0.0000 -0.0000 -0.7423 0.5369 1.0000 0.3650 0 0 0 0 0 0 0 0
0.0000 -0.0000 -0.0000 -0.0000 0.3712 -0.2685 -0.4619 0.0633 -0.0000 0.1232 -0.5000 1.0000 0.4513 0.8767 -0.0000 0.0000
0.0000 -0.0000 -0.0000 -0.0000 -0.6526 -0.7138 -0.3826 0.0722 -0.0000 0.1406 -0.4141 0.7002 0.3738 1.0000 -0.0000 -0.0000
0.0000 -0.0000 -0.0000 -0.0000 0.5000 0.2881 -0.0964 0.5121 -0.9876 -0.0000 0.0482 0.0000 0.5779 -0.0000 -0.0000 0.0000
0.0000 -0.0000 -0.0000 -0.0000 -0.4730 1.0000 -0.1582 1.0000 -1.0000 -0.0000 0.0791 0.0000 0.9481 -0.0000 -0.0000 0.0000
0.0000 -0.0000 -0.0000 -0.0000 -1.0000 -0.5762 -0.0648 -0.5314 0.0000 1.0000 -0.0702 0.0000 0.0633 0.4399 -0.0000 0.0000
0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000
0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000
0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0762 0.0000 -0.0000 0.0000 -1.0000 0.0000 0.9026 -0.0000 -0.0000 0.0000
0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 0.2458 -0.0000 0.1232 0.0000 1.0000 -0.0000 0.8767 0.0000 0.0000
zT(6,6) % the NMP zeros present also in G(s)
csi_zT = nT(1:14, 6); % zero input state directions
v_zT = nT (15:16, 6); % zero input directions
We can verify the presence of the zero we expected.
All this means that for every controller such that it is able to stabilize the plant G(s) and the resulting closed loop system exhibits the stability property, we will have a blocking NMP zero in 0.0128 In T(s).
Remembering that, in the case of no presence of disturbances or noises,
we can state that there will be a lack of control autority with respect to the references of the form
and the inizial states
, sinche the outputs with respect to such references will be identically equal to 0. Let's visualize it with a simulation. rT = [ones(length(t),1),ones(length(t),1)]; % inizializing the 2 references
rT(i,1) = exp(zi_nmp(2)*t(i))* v_zT(1,1); % imposing values and directions
rT(i,2) = exp(zi_nmp(2)*t(i))* v_zT(2,1);
y5 = lsim(T_sys,rT,t, csi_zT(1:14)); % outputs
plot(t, y5, 'LineWidth',2); grid
title("Limitation on T: blocking property of the zero z_{2} = 0.0128");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','NorthEast')
We can see that in such directions the reference can not be reproduced at steady state: the 2 outputs are always identically equal to zero (except for the problems of approximation that brings initially the outputs to reach values of the order of
). An additional limitation is derived from the fact that in the SISO
design we require: In a multivariable context, this is equivalent to require:
where
is the largest singular value of S(s). SInce
, we obtain an imposition of a bandwith limitation equal to :
in the case of z real, like this case . So, if we want to require a sensitivity peak to be not more than 2 as usual, and a very small error at steady state (A << 1), we obtain that the bandiwidth can not be greater than
. So in our design we have to take into account the fact that 
We can see from similations that only requiring a bandiwidth below this quantity the MATLAB function mixsyn() is able to find a stabilizing controller coherent with our specifications (
. If we try to require more bandiwith than this, we obtain a controller such that again the closed loop bandiwidth for S(s) is not greater than
, and consequently an higher γ, increasing as it increases the quantity of bandiwith required. (Remembering that γ is an index of how much the resulting controller is far from the point of view of the performances from our initial specifications).
om = logspace(-2,5,1000);
wp2 = (s/M2+ omb2)/ (s+omb2*A2);
wp22 = (s/M2+ omb2)/ (s+omb2*A2)* eye(2); % diagonal weight
[KHinf2, Ghinf2, gamma2 gopt2] = mixsyn(Gs_nmp, wp22, [], []);
gamma2 % we have a gamma very close to 1
Cloop2 = loopsens(Gs_nmp, KHinf2);
r_unit2 = [ones(length(t),1),ones(length(t),1)]; % dummy reference to check if the controller works
ycl2 = lsim(Cloop2.To, r_unit2, t);
plot(t, ycl2, 'LineWidth', 2), grid
title("System's steady state with simple H_{inf} controller (NMP)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
sigmaWS2 = sigma(wp22,om);
sigmaWS2Inv = sigma(1/wp22,om);
sigmaSo2 = sigma(Cloop2.So,om);
semilogx(om,1./sigmaWS2(2,:),om, sigmaWS2Inv(1,:),'--',... % plotting the minimum singular value of the weight , the max of the inverted weight, and the max of S(s)
om, sigmaSo2(1,:),'LineWidth',2), grid;
legend('1/s_{min}(wS2)','s_{max}(wS2^{-1})','s_{max}(S_o)','Location','northeast');
ylabel('dB','FontName','courier','FontSize',12)
xlabel('$\omega$','Interpreter','latex','FontSize',14)
We can see that requiring exactly
= 0.0064 we obtain a controller perfectly coherent with our specifications, such that the maximum singular value of the closed loop sensitivity function S(s) is perfectly upper bounded from the minimum singular value of the designed weight (or of the maximum singular value of the inverse of our weight). 4. Robustness Analysis
Minimum Phase case
Let us now consider some robustness issues, taking as a case of study a situation in which the parameters
are uncertain and vary in intervals between +10% and -10% of their nominal value. Initially we start to consider the system G(s) exhibiting Minimum Phase Properties. NMatlabRed = [0.8500 0.3250 0.0980];
NMatlabYellow = [0.929 0.694 0.125 ];
NMatlabBlue = [0 0.4470 0.7410];
NMatlabViolet = [0.4940 0.1840 0.5560];
NMatlabGreen = [0.4660 0.6740 0.1880];
NMatlabCyan = [0.3010 0.7450 0.9330];
NMatlabBordeaux = [0.6350 0.0780 0.1840];
k1_mp_u = ureal('K1',3.33,'Range',[(3.33-(3.33/10)), (3.33+(3.33/10))]); % cm^3/Vs
k2_mp_u = ureal('K2',3.35,'Range',[(3.35-(3.35/10)), (3.35+(3.35/10))]); % cm^3/Vs
g1_mp_u = ureal('G1',0.7,'Range',[(0.7-(0.7/10)), (0.7+(0.7/10))]);
g2_mp_u = ureal('G2',0.6,'Range',[(0.6-(0.6/10)), (0.6+(0.6/10))]);
Gs_mp_unc = [(g1_mp_u*k1_mp_u*n1_mp/(1+t1_mp*s)), ((1-g2_mp_u)*k2_mp_u*n1_mp)/ ((1+t3_mp*s)*(1+t1_mp*s)); ((1-g1_mp_u)*k1_mp_u*n2_mp)/ ((1+t2_mp*s)*(1+t4_mp*s)), (g2_mp_u*k2_mp_u*n2_mp)/(1+t2_mp*s)]
Uncertain continuous-time state-space model with 2 outputs, 2 inputs, 6 states.
The model uncertainty consists of the following blocks:
G1: Uncertain real, nominal = 0.7, range = [0.63,0.77], 2 occurrences
G2: Uncertain real, nominal = 0.6, range = [0.54,0.66], 2 occurrences
K1: Uncertain real, nominal = 3.33, range = [3,3.66], 2 occurrences
K2: Uncertain real, nominal = 3.35, range = [3.02,3.69], 2 occurrences
Model Properties
Type "Gs_mp_unc.NominalValue" to see the nominal value and "Gs_mp_unc.Uncertainty" to interact with the uncertain elements.
We want to represent this parametric uncertainty as an unstructured input multiplicative uncertainty of the form:
with
and thw weight
such that it satisfy the requirement that:
, with
the set of all possible perturbed plants.We compute a 3° order diagonal weight through the command ucover().
w =logspace (-3,2, 1000); % frequency span
Gs_mp_unc_sampled = gridureal(Gs_mp_unc,30);
G_star_mp = inv(Gs_mp)* (Gs_mp_unc_sampled- Gs_mp); % showing the singular values of the matrix G^-1( G_p - G)
[PwM1,InfoM1] = ucover(Gs_mp_unc_sampled,Gs_mp,3, 'InputMult'); % computing a suitable weight
[m, p] = bode(InfoM1.W1, w);
plot_weight = semilogx(w, 20*log10(m(1,:)));
set(plot_weight,'LineWidth',4, 'Color',NMatlabRed), grid on;
Leg2 = legend('$\sigma$($G^{-1}$($G_{u}$-G)','$w_m$ with n = 3','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
title('Weight for singular values, multiplicative uncertanity','FontSize',10);
Let us visualize the resulting input multiplicative representation of the uncertainty in terms of relative error:
[PwM1,InfoM1] = ucover(Gs_mp_unc_sampled,Gs_mp,[3 3], 'InputMult'); % choiche of a diagonal weight
Rel_P = (Gs_mp- Gs_mp_unc_sampled )/Gs_mp; % relative error
bodemag(Rel_P,'b--',InfoM1.W1,'r'), grid on;
Leg2 = legend('sampled','$w_m$ with n = 3','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
title('Relative error, multiplicative uncertanity','FontSize',10);
bodemag(PwM1,'b',Gs_mp,'r'), grid on;
Leg2 = legend('uncertain','nominal','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
title('multiplicative uncertanity representation','FontSize',10);
We can see from this plots that the most uncertain situation is at low frequency, while at high frequency the diagonal elements of the uncertainty representations are both very close to the nominal process behavior .This result is also highlighted in the Vadigepalli, Gatzke, Doyle's paper.
Let us now compute two different stabilizing decentralized controllers for the nominal case. We want to compare the performance of two different controllers in terms of reference reproduction's capabilities and also in terms of robust stability properties.
The first controller is the decentralized controller PI proposed in the Johansson's paper, the second one is another diagonal controller designed with an
norm control approach.
% controller for the G(s) with MP property
PI1_mp = K1_mp*(1+ (1/(T1_mp*s)));
PI2_mp = K2_mp*(1+ (1/(T2_mp*s)));
PI_mp = [ PI1_mp ,0; 0, PI2_mp] % decentralized PI controller with Johansson's parameters
PI_mp =
From input 1 to output...
90 s + 3
1: --------
30 s
2: 0
From input 2 to output...
1: 0
108 s + 2.7
2: -----------
40 s
Continuous-time transfer function.
Model Properties
Cloop3 = loopsens(Gs_mp, PI_mp);
r_unit2 = [ones(length(t),1),ones(length(t),1)];
ycl3 = lsim(Cloop3.To, r_unit2, t);
plot(t, ycl3, 'LineWidth', 2), grid
title("System's steady state with Decentralized PI controller (MP)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
On the other hand, we need for comparison an
controller for the MP case. We will use the input-output paring derived to be better in the Relative Array Gain analysis:
. This implies that the performance weight of the decentralized
controller will have off-diagonal components
and the diagonal ones = 0.Since we don't have costraints on the control effort, we just design a weight for the sensitivity function S(s).
om = logspace(-2,5,1000);
wp3 = (s/M3+ omb3)/ (s+omb3*A3);
wp33 = (s/M3+ omb3)/ (s+omb3*A3)* [0,1;1,0]; % off-diagonal weight
[KHinf3, Ghinf3, gamma3 gopt2] = mixsyn(Gs_mp, wp33, [], []);
gamma3 % we have a gamma very close to 1
Cloop5 = loopsens(Gs_mp, KHinf3);
ycl5 = lsim(Cloop5.To, r_unit2, t);
plot(t, ycl5, 'LineWidth', 2), grid
title("System's steady state with H_{inf} controller (MP)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
We can see the differences in terms of nominal performance: the
controller allows a faster transient and not overshoot.
Let us check if this 2 controllers are robust stable, i.e. if they satisfy the necessary and sufficient condition
.
L_MP_PI = PI_mp*PwM1; % open loop function PI
L_MP_Hinf = KHinf3*PwM1; % open loop function H-infinity
CL_MP_PI = feedback(L_MP_PI, eye(2)); % closed loop function PI
CL_MP_Hinf= feedback(L_MP_Hinf, eye(2)); % closed loop function H-infinity
[stabmarg1,wcu1] = robstab(CL_MP_PI)
stabmarg1 =
LowerBound: 1.7088
UpperBound: 1.7123
CriticalFrequency: 1.0000e-04
wcu1 =
Gs_mp_unc_sampled_InputMultDelta: [2×2 ss]
[stabmarg2,wcu2] = robstab(CL_MP_Hinf)
stabmarg2 =
LowerBound: 1.7102
UpperBound: 1.7137
CriticalFrequency: 1.0000e-04
wcu2 =
Gs_mp_unc_sampled_InputMultDelta: [2×2 ss]
Being the upper and lower bounds striclty >1 in both case, we can state that both controllers are able to robust stabilize the uncertain system with MP property.
Let us check if they can also guarantee robust performance, i.d. if 
where M is the
component of the augmented plant
the only possible source of instability of the overall control structured system
and
is the structured singular value.
[perfmarg_MP_Hinf,wcu,report,info1] = robustperf(CL_MP_Hinf)
perfmarg_MP_Hinf =
LowerBound: 0.9642
UpperBound: 0.9643
CriticalFrequency: 0.0310
wcu =
Gs_mp_unc_sampled_InputMultDelta: [2×2 ss]
report =
'Uncertain system does not achieve performance robustness to modeled uncertainty. '
' -- The tradeoff of model uncertainty and system gain is balanced at a level of 96.4% of the modeled uncertainty. '
' -- A model uncertainty of 96.4% can lead to input/output gain of 1.04 at 0.031 rad/seconds. '
' -- Sensitivity with respect to each uncertain element is: '
' 13% for Gs_mp_unc_sampled_InputMultDelta. Increasing Gs_mp_unc_sampled_InputMultDelta by 25% decreases the margin by 3.25%.'
info1 =
Sensitivity: [1×1 struct]
Frequency: [113×1 double]
BadUncertainValues: [113×1 struct]
MussvBnds: [1×2 frd]
MussvInfo: [1×1 struct]
[perfmarg_MP_PI,wcu,report,info2] = robustperf(CL_MP_PI)
perfmarg_MP_PI =
LowerBound: 0.8010
UpperBound: 0.8022
CriticalFrequency: 0.0567
wcu =
Gs_mp_unc_sampled_InputMultDelta: [2×2 ss]
report =
'Uncertain system does not achieve performance robustness to modeled uncertainty. '
' -- The tradeoff of model uncertainty and system gain is balanced at a level of 80.1% of the modeled uncertainty. '
' -- A model uncertainty of 80.2% can lead to input/output gain of 1.25 at 0.0567 rad/seconds. '
' -- Sensitivity with respect to each uncertain element is: '
' 15% for Gs_mp_unc_sampled_InputMultDelta. Increasing Gs_mp_unc_sampled_InputMultDelta by 25% decreases the margin by 3.75%.'
info2 =
Sensitivity: [1×1 struct]
Frequency: [69×1 double]
BadUncertainValues: [69×1 struct]
MussvBnds: [1×2 frd]
MussvInfo: [1×1 struct]
While we had robust stability for both control systems, we don't have robust performance capabilities, being the lower and the upper bounds for robust performance strictly <1 in both cases.
Let us visualize how the structured singulare value varies with variation of frequencies.
% extracting the values for the H inf controller
a = info1.MussvBnds(1,1);
b = info1.MussvBnds(1,2);
max_mp_Hinf = zeros(length(info1.Frequency),1);
for i=1:size(info1.Frequency)
max_mp_Hinf(i)= a.Response(:,:,i);
min_mp_Hinf = zeros(length(info1.Frequency),1);
for i=1:length(info1.Frequency)
min_mp_Hinf(i)= b.Response(:,:,i);
% extracting the values for the PI controller
c = info2.MussvBnds(1,1);
d = info2.MussvBnds(1,2);
max_mp_PI = zeros(length(info2.Frequency),1);
for i=1:length(info2.Frequency)
max_mp_PI(i)= c.Response(:,:,i);
min_mp_PI = zeros(length(info2.Frequency),1);
for i=1:length(info2.Frequency)
min_mp_PI(i)= d.Response(:,:,i);
semilogx(info1.Frequency, max_mp_Hinf,'b--',info1.Frequency, min_mp_Hinf,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
title('µ bounds H-inf CL');
semilogx(info2.Frequency, max_mp_PI,'b--',info2.Frequency, min_mp_PI,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
This plot confirms our expectations: both controllers can not achieve robust performance since the structured singular value is >1 for some ω, but we can see that the
one is significantly closer to achieve RP: the upper and lower bounds are higher than the ones of the PI closed loop and also the condition
< 1 is vialoted at very few frequencies.
Let us try to design a
controller able to guarantee robust performance trhough μ-synthesis control. We have to build the augumented plant and run the D-K algorithm.
wp44 = ((s/M4+ omb4)/ (s+omb4*A4))* [0,1;1,0]; % off diagonal weight
% inizializing the inputs and the outputs of each block
PwM1.u = 'u'; % control input
PwM1.y = 'y'; % not perturbed output
wp44.y = 'z'; % performance variable
Sum3 = sumblk('v = r - y', 2); % measured variables
Pe_MP = connect(PwM1,wp44,Sum3, {'r','u'},{'z','v' }); % augmented plant
ncont = 2; % two control inputs
nmeas = 2; % two measurements
[Krob1,rpMU1] = musyn(Pe_MP ,nmeas,ncont)
D-K ITERATION SUMMARY:
-----------------------------------------------------------------
Robust performance Fit order
-----------------------------------------------------------------
Iter K Step Peak MU D Fit D
1 30.04 30.04 30.28 8
2 10.88 10.88 10.88 8
3 3.278 2.964 2.983 20
4 2.784 1.797 1.807 12
5 1.704 1.454 1.46 12
6 1.422 1.301 1.306 12
7 1.28 1.204 1.212 20
8 1.2 1.144 1.153 20
9 1.135 1.094 1.103 20
10 1.08 1.05 1.059 20
Best achieved robust performance: 1.05
Krob1 =
A =
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 x27 x28 x29 x30 x31 x32 x33 x34
x1 -0.01328 1 6.232e-10 2.358e-10 1.329e-10 3.849e-14 -1.685e-11 -6.507e-10 -2.576e-10 -1.384e-10 -1.029e-09 2.088e-10 1.229e-09 -2.482e-10 -9.192e-10 1.019e-09 -1.712e-10 -6.92e-11 -2.254e-10 1.551e-10 6.101e-11 2.05e-10 1.15e-08 -1.446e-08 -1.529e-14 6.513e-12 2.574e-10 9.181e-11 5.624e-11 1.411e-14 -6.451e-12 -2.522e-10 -9.273e-11 -5.493e-11
x2 -0.000435 0.1163 7.053 3.388 1.283 0.0003546 -0.08741 -3.66 -1.222 -0.8192 -10.87 2.009 6.158 -2.254 -9.522 4.464 -1.551 -0.6273 -2.051 1.507 0.5947 1.99 196.5 -0.09352 -0.0001713 0.06759 2.63 1.096 0.5263 0.0001599 -0.04579 -1.889 -0.6306 -0.4227
x3 -0.01262 4.378 171.7 114.5 43.34 0.01198 -2.953 -123.6 -41.29 -27.68 -367.3 67.87 208 -76.16 -321.7 150.8 -52.41 -21.19 -69.28 50.92 20.09 67.24 6640 -3.159 -0.005786 2.284 88.86 37.03 17.78 0.005402 -1.547 -63.8 -21.31 -14.28
x4 -3.859e-18 9.36e-16 3.563e-14 -104.5 86.75 4.058e-18 1.101e-15 3.026e-14 2.207e-14 4.621e-15 -1.871e-13 2.507e-14 -9.62e-14 -2.316e-14 -1.559e-13 -1.075e-13 -1.54e-14 -6.227e-15 -2.059e-14 1.977e-14 7.848e-15 2.604e-14 3.636e-12 4.29e-12 -2.484e-18 4.505e-16 1.888e-14 8.151e-15 3.787e-15 2.709e-18 -1.044e-16 -8.243e-15 -2.272e-16 -2.299e-15
x5 -0.01034 3.589 140.8 -37.4 -64.38 0.00982 -2.42 -101.3 -33.84 -22.68 -301.1 55.63 170.5 -62.43 -263.7 123.6 -42.96 -17.37 -56.79 41.74 16.47 55.11 5442 -2.59 -0.004742 1.872 72.83 30.36 14.57 0.004428 -1.268 -52.3 -17.46 -11.71
x6 2.58e-13 -1.132e-10 -4.484e-09 -1.532e-09 -9.953e-10 -0.01328 1 6.224e-09 2.694e-09 1.282e-09 4.474e-09 -1.203e-09 -1.265e-08 1.581e-09 4.252e-09 -1.112e-08 1.103e-09 4.46e-10 1.445e-09 -8.668e-10 -3.394e-10 -1.148e-09 -1.468e-08 2.209e-07 7.324e-14 -4.507e-11 -1.757e-09 -5.942e-10 -3.916e-10 -5.808e-14 5.246e-11 1.956e-09 7.832e-10 4.145e-10
x7 0.0002309 -0.09807 -3.906 -1.304 -0.8743 -0.0002789 0.1444 7.787 3.75 1.494 3.458 -1.1 -11.67 1.429 3.392 -10.27 0.9931 0.4015 1.301 -0.7984 -0.3126 -1.057 -0.08686 241.9 6.595e-05 -0.03976 -1.554 -0.519 -0.3479 -5.961e-05 0.06113 2.235 0.989 0.4477
x8 0.007799 -3.313 -132 -44.06 -29.54 -0.007343 5.326 196.5 126.7 50.47 116.8 -37.15 -394.2 48.29 114.6 -346.9 33.55 13.56 43.96 -26.97 -10.56 -35.72 -2.934 8172 0.002228 -1.343 -52.51 -17.54 -11.75 -0.002014 2.065 75.51 33.41 15.12
x9 5.032e-18 -2.445e-15 -9.731e-14 -3.06e-14 -2.223e-14 -4.398e-18 4.421e-15 1.601e-13 -104.5 86.75 4.948e-14 -2.123e-14 -3.363e-13 3.098e-14 5.391e-14 -3.029e-13 2.186e-14 8.836e-15 2.851e-14 -1.476e-14 -5.747e-15 -1.96e-14 7.871e-13 6.695e-12 9.919e-19 -9.492e-16 -3.659e-14 -1.182e-14 -8.286e-15 -5.543e-19 1.237e-15 4.475e-14 1.888e-14 9.308e-15
x10 0.006393 -2.716 -108.2 -36.12 -24.21 -0.006019 4.366 161.1 -27.37 -58.54 95.77 -30.45 -323.1 39.58 93.94 -284.3 27.5 11.12 36.03 -22.11 -8.657 -29.28 -2.405 6699 0.001826 -1.101 -43.04 -14.37 -9.634 -0.001651 1.693 61.89 27.39 12.4
x11 -0.02447 8.644 338.4 118.9 96.35 0.02375 -5.985 -250 -83.47 -55.95 -756.6 133.6 412.3 -148.6 -658.5 299.5 -101.6 -41.07 -134.2 101 39.82 133.3 1.373e+04 -6.46 -0.01077 3.657 146.2 50.06 37.07 0.01058 -3.098 -127.5 -42.58 -28.55
x12 -0.02446 8.644 338.4 118.9 96.35 0.02375 -5.985 -250 -83.47 -55.95 -756.6 133.5 412.3 -148.6 -658.4 299.5 -101.6 -41.06 -134.2 100.9 39.82 133.3 1.373e+04 -6.678 -0.01077 3.657 146.2 50.06 37.07 0.01058 -3.098 -127.5 -42.58 -28.54
x13 -8.713e-07 0.0002955 0.01154 0.004915 0.002331 8.55e-07 -0.0001538 -0.00688 -0.00195 -0.001603 -0.02889 0.0205 0.01 -0.005283 -0.02494 0.006176 -0.003599 -0.001455 -0.004764 0.003707 0.001464 0.004894 0.4412 0.1244 -4.065e-07 0.0001266 0.005086 0.001922 0.001085 4.071e-07 -9.919e-05 -0.004197 -0.001328 -0.0009527
x14 0.01472 -6.638 -263.3 -87.92 -58.94 -0.01324 10.48 384.8 142.8 109.7 225.9 -66.9 -791.2 90.22 220 -701.3 63.19 25.54 82.68 -47.79 -18.69 -63.34 -5.998 1.68e+04 0.00384 -2.625 -102 -34.06 -22.83 -0.002866 3.157 116.7 41.98 31.08
x15 9.9e-07 -0.0004812 -0.01915 -0.00602 -0.004376 -8.652e-07 0.0008701 0.03152 0.0142 0.006394 0.009721 -0.004175 -0.0662 0.03735 0.0106 -0.05963 0.004301 0.001738 0.005609 -0.002904 -0.00113 -0.003855 0.1553 1.318 1.95e-07 -0.0001868 -0.0072 -0.002325 -0.00163 -1.088e-07 0.0002434 0.008806 0.003716 0.001832
x16 0.01472 -6.639 -263.3 -87.93 -58.94 -0.01324 10.48 384.9 142.8 109.7 225.9 -66.91 -791.3 90.29 220 -701.5 63.19 25.55 82.69 -47.8 -18.69 -63.34 -5.783 1.681e+04 0.003841 -2.625 -102 -34.07 -22.83 -0.002867 3.157 116.7 41.99 31.09
x17 0.000617 0.1222 3.429 1.765 0.9 0.0005514 -0.4215 -16.32 -5.451 -3.652 -79.14 0.9272 9.275 1.744 -60.04 8.126 2.664 1.081 3.619 2.346 0.9055 3.121 1543 -1.011 0.001259 -1.8 -61.53 -49.43 4.797 -3.235e-05 -0.1402 -5.188 -1.733 -1.161
x18 -0.0003107 -0.06033 -1.685 -0.87 -0.4439 -0.0002735 0.2103 8.14 2.72 1.822 39.49 -0.4511 -4.611 -0.8856 29.96 -4.046 -1.328 -0.5698 -1.822 -1.163 -0.449 -1.548 -770.2 0.8034 -0.0006303 0.9004 30.78 24.72 -2.396 1.722e-05 0.06986 2.583 0.8633 0.5782
x19 7.805e-05 0.01496 0.4162 0.2154 0.11 6.802e-05 -0.0525 -2.031 -0.6789 -0.4546 -9.862 0.1108 1.147 0.2236 -7.48 1.007 0.3336 0.1494 0.4321 0.2893 0.1117 0.385 192.4 -0.2097 0.0001577 -0.2251 -7.696 -6.181 0.5985 -4.465e-06 -0.01742 -0.644 -0.2152 -0.1441
x20 -0.001277 -0.3061 -9.79 -3.271 -2.192 0.002551 0.07231 -1.456 1.614 -0.4688 -7.173 13.09 -39.49 -9.731 -10.49 -45.98 -5.745 -2.325 -7.762 11.15 4.421 14.72 -0.3622 1697 -0.00116 0.02189 2.183 0.7279 0.4886 0.002374 -1.959 -68.9 -51.61 3.107
x21 0.0006396 0.1525 4.872 1.629 1.091 -0.001276 -0.03503 0.7682 -0.7886 0.2425 3.592 -6.549 19.66 4.873 5.253 22.91 2.877 1.164 3.887 -5.566 -2.237 -7.363 0.4968 -846.8 0.0005804 -0.01116 -1.1 -0.3666 -0.2462 -0.001187 0.9799 34.46 25.81 -1.551
x22 -0.0001599 -0.03809 -1.217 -0.4069 -0.2723 0.0003191 0.008674 -0.195 0.1957 -0.06121 -0.8968 1.637 -4.908 -1.218 -1.312 -5.722 -0.7195 -0.2912 -0.9721 1.392 0.5655 1.816 -0.1778 211.5 -0.0001451 0.0028 0.2754 0.09175 0.06164 0.0002968 -0.245 -8.615 -6.452 0.3876
x23 -8.047e-16 3.318e-13 1.31e-11 4.716e-12 2.852e-12 7.467e-16 -4.237e-13 -1.595e-11 -6.654e-12 -3.332e-12 -1.732e-11 3.95e-12 3.144e-11 -4.917e-12 -1.585e-11 2.701e-11 -3.41e-12 -1.379e-12 -4.479e-12 2.894e-12 1.136e-12 3.829e-12 -0.4 -4.761e-10 -2.674e-16 1.343e-13 5.273e-12 1.833e-12 1.163e-12 2.327e-16 -1.445e-13 -5.512e-12 -2.121e-12 -1.183e-12
x24 2.499e-16 -1.029e-13 -4.063e-12 -1.464e-12 -8.84e-13 -2.319e-16 1.309e-13 4.932e-12 2.055e-12 1.03e-12 5.4e-12 -1.228e-12 -9.714e-12 1.527e-12 4.937e-12 -8.339e-12 1.059e-12 4.279e-13 1.391e-12 -8.999e-13 -3.533e-13 -1.191e-12 -4.455e-11 -0.4 8.326e-17 -4.166e-14 -1.636e-12 -5.689e-13 -3.608e-13 -7.258e-17 4.476e-14 1.708e-12 6.565e-13 3.667e-13
x25 -0.00146 0.4749 18.5 8.174 3.665 0.001447 -0.1614 -8.245 -1.608 -2.052 -51.59 8.352 9.167 -8.838 -44.23 3.064 -6 -2.426 -7.953 6.392 2.526 8.435 829.9 410.4 -0.014 1.206 8.323 3.201 1.763 0.0007309 -0.1481 -6.473 -1.92 -1.492
x26 -0.0001797 0.07268 2.796 1.188 0.5746 0.0002135 -0.0533 -2.256 -0.6986 -0.5149 -9.072 1.14 2.932 -1.156 -7.575 2.003 -0.7393 -0.2986 -0.9762 0.9192 0.3622 1.214 147.6 19.61 -0.0001121 -0.0428 1.134 0.4325 0.2431 9.017e-05 -0.02814 -1.157 -0.3748 -0.2611
x27 0.0002828 0.2757 9.298 4.213 2.127 0.0009954 -0.5629 -22.09 -7.38 -4.944 -99.66 3.198 16.97 -0.4936 -76.99 13.86 1.308 0.5409 1.796 4.233 1.648 5.616 1879 -1.572 0.001214 -1.928 -65.6 0.5666 0.2325 0.0001365 -0.2053 -7.805 -2.608 -1.747
x28 -0.0009465 0.3345 13.1 5.38 2.692 0.0009189 -0.2316 -9.675 -3.231 -2.165 -29.23 5.167 15.96 -5.748 -25.44 11.6 -3.929 -1.589 -5.193 3.906 1.541 5.158 418.3 -0.4068 -0.0004168 0.1415 5.656 -102.4 87.96 0.0004093 -0.1199 -4.934 -1.648 -1.105
x29 0.003771 -1.025 -41.35 -16.66 -8.324 -0.00262 0.4047 18.07 6.032 4.045 27.63 -16.7 -45.76 21.09 32.04 -31.99 15.77 6.384 20.89 -11.14 -4.41 -14.69 -24.05 -0.04842 0.002554 -2.109 -74.88 -138.6 -104.3 -0.001419 0.28 12.05 4.024 2.698
x30 0.001913 -0.8858 -35.17 -11.5 -7.932 -0.001704 1.471 53.72 23.8 10.97 25.8 -8.483 -111.3 11.75 25.86 -99.26 8.248 3.334 10.78 -6.01 -2.346 -7.968 100.7 2104 0.0004579 -0.3479 -13.48 -4.451 -3.03 -0.0136 1.431 15.79 6.512 3.312
x31 5.275e-05 -0.05776 -2.224 -0.6867 -0.5111 2.08e-06 0.09425 3.262 1.56 0.6531 0.3984 0.04637 -8.304 0.265 0.442 -7.861 0.2255 0.09104 0.2842 0.1043 0.04374 0.1342 23.18 199.8 -2.269e-05 -0.01765 -0.6337 -0.1999 -0.1447 1.007e-05 -0.05235 0.6119 0.3038 0.1207
x32 -0.00119 -0.4359 -14.7 -4.912 -3.29 0.002617 0.2347 4.087 4.276 0.6501 -4.548 13.44 -55.24 -9.395 -8.297 -61.07 -5.4 -2.186 -7.336 11.61 4.614 15.28 -0.9056 2085 -0.001223 -0.01475 0.8966 0.2979 0.2006 0.002576 -2.114 -74.33 -1.965 -1.785
x33 0.0004843 -0.2185 -8.665 -2.894 -1.94 -0.0004356 0.3448 12.66 5.553 2.596 7.44 -2.201 -26.01 2.971 7.245 -23.05 2.079 0.8405 2.721 -1.573 -0.6148 -2.084 -0.359 475.1 0.0001264 -0.08637 -3.357 -1.121 -0.7513 -9.435e-05 0.1039 3.838 -102.9 87.56
x34 -0.002861 0.4932 21.69 7.239 4.854 0.00384 -1.15 -45.96 -18.12 -9.577 -32.66 19.58 56 -19.26 -34.98 39.72 -12.52 -5.063 -16.6 15.64 6.174 20.64 0.01414 -141.6 -0.001494 0.3241 13.8 4.608 3.089 0.002478 -2.137 -75.82 -138.9 -104.5
B =
u1 u2
x1 -9.5e-10 9.238e-10
x2 -4.467e-11 2.183e-11
x3 -4.602e-12 -1.894e-12
x4 -3.006e-13 -2.744e-13
x5 -1.354e-13 4.559e-14
x6 1.205e-09 -1.413e-08
x7 -2.223e-11 -2.847e-11
x8 -1.218e-13 1.665e-12
x9 -6.529e-14 -4.28e-13
x10 -8.181e-13 -5.564e-13
x11 -99.36 0.002303
x12 -99.29 0.01625
x13 -0.03647 -0.007969
x14 0.005347 -88.57
x15 -0.01288 -0.08427
x16 -0.01248 -88.71
x17 -217.4 0.03289
x18 108.6 -0.03552
x19 -27.14 0.009449
x20 -0.00229 -194.2
x21 -0.02502 96.96
x22 0.01069 -24.23
x23 8 3.044e-11
x24 3.677e-12 8
x25 -68.61 -26.27
x26 -15.14 -1.258
x27 -254.5 0.05914
x28 -34.57 0.01106
x29 -79.35 0.02508
x30 -8.399 -134.5
x31 -1.922 -15.57
x32 0.03168 -227.8
x33 0.01269 -30.38
x34 0.02958 -68.38
C =
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 x27 x28 x29 x30 x31 x32 x33 x34
y1 -0.09898 34.35 1347 472.3 383.8 0.09399 -23.17 -970 -323.9 -217.1 -2882 532.5 1632 -597.6 -2524 1183 -411.2 -166.3 -543.6 399.5 157.6 527.5 5.209e+04 -24.79 -0.04539 17.92 697.2 290.6 139.5 0.04238 -12.14 -500.6 -167.2 -112.1
y2 0.06119 -25.99 -1035 -345.7 -231.8 -0.05761 41.79 1542 568.4 439.7 916.7 -291.5 -3093 378.9 899.2 -2722 263.2 106.4 344.9 -211.6 -82.86 -280.3 -23.02 6.412e+04 0.01748 -10.54 -412 -137.6 -92.22 -0.0158 16.2 592.4 262.1 118.7
D =
u1 u2
y1 0 0
y2 0 0
Continuous-time state-space model.
Model Properties
rpMU1 = 1.0502
After 10 iterations, the best attempt of musyn() using D-K algorithm ends with a
order controller able to guarantee a maximum structured singular value of 1.052, so not robust performant in absolute terms, but very close. We expect lower and upper bounds for robust performance indices very close to one.
L_MP_rob = Krob1*PwM1; % open loop function robust performance controller (MP)
CL_MP_rob = feedback(L_MP_rob, eye(2)); % closed loop robust performance controller (MP)
[perfmarg_MP_rob,wcu,report,info3] = robustperf(CL_MP_rob)
perfmarg_MP_rob =
LowerBound: 0.9805
UpperBound: 0.9806
CriticalFrequency: 12.3801
wcu =
Gs_mp_unc_sampled_InputMultDelta: [2×2 ss]
report =
'Uncertain system does not achieve performance robustness to modeled uncertainty. '
' -- The tradeoff of model uncertainty and system gain is balanced at a level of 98.1% of the modeled uncertainty. '
' -- A model uncertainty of 98.1% can lead to input/output gain of 1.02 at 12.4 rad/seconds. '
' -- Sensitivity with respect to each uncertain element is: '
' 14% for Gs_mp_unc_sampled_InputMultDelta. Increasing Gs_mp_unc_sampled_InputMultDelta by 25% decreases the margin by 3.5%.'
info3 =
Sensitivity: [1×1 struct]
Frequency: [121×1 double]
BadUncertainValues: [121×1 struct]
MussvBnds: [1×2 frd]
MussvInfo: [1×1 struct]
Let us visualize the behavior of the structured singular value in a suitable range of frequencies.
% extracting the structured singular values for the µ-syn robust controller
p = info3.MussvBnds(1,1);
o = info3.MussvBnds(1,2);
max_mp_rob = zeros(length(info3.Frequency),1);
for i=1:length(info3.Frequency)
max_mp_rob(i)= p.Response(:,:,i);
min_mp_rob = zeros(length(info3.Frequency),1);
for i=1:length(info3.Frequency)
min_mp_rob(i)= o.Response(:,:,i);
semilogx(info3.Frequency, max_mp_rob,'b--',info3.Frequency, min_mp_rob,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
title('µ bounds robust performant CL');
We can see that µ is always <1 exept for a small neighborood of the critical frequency
. This result is reasonably the best achieved in terms of robust performance, better than the
nominal controller which achieves a robust performance index around 96.42%, and than the PI that achieves about 80%. Let us compute the step response and try to visualize some possible behaviors (20) of closed loop system under the µ-syn controller.
L_MP_rob = Krob1* PwM1; % robust open loop
CL_MP_rob = feedback(L_MP_rob, eye(2)); % robust closed loop
samples_MP = usample(CL_MP_rob ,20);
lsim(samples_MP , 'b',r_unit2, t);
title("System's steady state with µ-synthesis controller (MP)");
Non-Minimum Phase case
k1_nmp_u = ureal('K1',3.14,'Range',[(3.14-(3.14/10)), (3.14+(3.14/10))]); % cm^3/Vs
k2_nmp_u = ureal('K2',3.29,'Range',[(3.29-(3.29/10)), (3.29+(3.29/10))]); % cm^3/Vs
g1_nmp_u = ureal('G1',0.43,'Range',[(0.43-(0.43/10)), (0.43+(0.43/10))]);
g2_nmp_u = ureal('G2',0.34,'Range',[(0.34-(0.34/10)), (0.34+(0.34/10))]);
Gs_nmp_unc = [(g1_nmp_u*k1_nmp_u*n1_nmp/(1+t1_nmp*s)), ((1-g2_nmp_u)*k2_nmp_u*n1_nmp)/ ((1+t3_nmp*s)*(1+t1_nmp*s)); ((1-g1_nmp_u)*k1_nmp_u*n2_nmp)/ ((1+t2_nmp*s)*(1+t4_nmp*s)), (g2_nmp_u*k2_nmp_u*n2_nmp)/(1+t2_nmp*s)]
Uncertain continuous-time state-space model with 2 outputs, 2 inputs, 6 states.
The model uncertainty consists of the following blocks:
G1: Uncertain real, nominal = 0.43, range = [0.387,0.473], 2 occurrences
G2: Uncertain real, nominal = 0.34, range = [0.306,0.374], 2 occurrences
K1: Uncertain real, nominal = 3.14, range = [2.83,3.45], 2 occurrences
K2: Uncertain real, nominal = 3.29, range = [2.96,3.62], 2 occurrences
Model Properties
Type "Gs_nmp_unc.NominalValue" to see the nominal value and "Gs_nmp_unc.Uncertainty" to interact with the uncertain elements.
We want to represent this parametric uncertainty as an unstructured input multiplicative uncertainty of the form:
with
and thw weight
such that it satisfy the requirement that:
, with
the set of all possible perturbed plants.We compute a 3° order diagonal weight through the command ucover().
w =logspace (-3,2, 1000); % frequency span
Gs_nmp_unc_sampled = gridureal(Gs_nmp_unc,30);
G_star_nmp = inv(Gs_nmp)* (Gs_nmp_unc_sampled- Gs_nmp); % showing the singular values of the matrix G^-1( G_p - G)
[PwM2,InfoM2] = ucover(Gs_nmp_unc_sampled,Gs_nmp,3, 'InputMult'); % computing a suitable weight
[m, p] = bode(InfoM2.W1, w);
plot_weight = semilogx(w, 20*log10(m(1,:)));
set(plot_weight,'LineWidth',4, 'Color',NMatlabRed), grid on;
Leg2 = legend('$\sigma$($G^{-1}$($G_{u}$-G)','$w_m$ with n = 3','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
title('Weight for singular values, multiplicative uncertanity','FontSize',10);
Let us visualize the resulting input multiplicative representation of the uncertainty in terms of relative error:
[PwM2,InfoM2] = ucover(Gs_nmp_unc_sampled,Gs_nmp,[3 3], 'InputMult'); % choiche of a diagonal weight
Rel_P = (Gs_nmp- Gs_nmp_unc_sampled )/Gs_nmp; % relative error
bodemag(Rel_P,'b--',InfoM2.W1,'r'), grid on;
Leg2 = legend('sampled','$w_m$ with n = 3','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
title('Relative error, multiplicative uncertanity','FontSize',10);
bodemag(PwM2,'b',Gs_nmp,'r'), grid on;
Leg2 = legend('uncertain','nominal','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
title('multiplicative uncertanity representation','FontSize',10);
We can see from this plots that the most uncertain situation is at low frequency, while at high frequency the diagonal elements of the uncertainty representations are both very close to the nominal process behavior .This result is also highlighted in the Vadigepalli, Gatzke, Doyle's paper.
Let us now compute two different stabilizing decentralized controllers for the nominal case.
The first controller is the decentralized controller PI proposed in the Johansson's paper, the second one is another diagonal controller designed with an
norm control approach.
% controller for the G(s) with NMP property
PI1_nmp = K1_nmp*(1+ (1/(T1_nmp*s)));
PI2_nmp = K2_nmp*(1+ (1/(T2_nmp*s)));
PI_nmp = [PI1_nmp, 0; 0,PI2_nmp] % decentralized PI controller with Johansson's parameters
PI_nmp =
From input 1 to output...
165 s + 1.5
1: -----------
110 s
2: 0
From input 2 to output...
1: 0
-26.4 s - 0.12
2: --------------
220 s
Continuous-time transfer function.
Model Properties
Cloop4 = loopsens(Gs_nmp, PI_nmp);
t_nmp = 1:0.001:4000; % inizializing the time span for the non minimum phase transient!
r_unit3 = [ones(length(t_nmp),1),ones(length(t_nmp),1)];
ycl4 = lsim(Cloop4.To, r_unit3, t_nmp);
plot(t_nmp, ycl4, 'LineWidth', 2), grid
title("System's steady state with Decentralized PI controller (NMP)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
Note the difference in terms of time scale: the PI controller produces an output of the system with NMP property with a transient about 10 times slower than in the MP case.
For the
controller we can use the one designed in the previous section for the NMP case, the one that takes into account the bandwith limitations imposed by the presence of the
transmittion zero.
om = logspace(-2,5,1000);
wp55 = (s/M7+ omb7)/ (s+omb7*A7)* [0,1;1,0]; % off-diagonal weight
[KHinf7, Ghinf7, gamma7, gopt7] = mixsyn(Gs_nmp, wp55, [], []);
gamma7 % we have a gamma very close to 1
Cloop7 = loopsens(Gs_nmp, KHinf7);
ycl7 = lsim(Cloop7.To, r_unit2, t);
plot(t, ycl7, 'LineWidth', 2), grid
title("System's steady state with H_{inf} controller (MP)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
We can see the differences in terms of nominal performance: the
controller allows a faster transient and not overshoot. In both situations the outputs are slower than the corresponding MP case of about a order of magnitude.
Let us check if this 2 controllers are robust stable, i.e. if they satisfy the necessary and sufficient condition
.
L_NMP_PI = PI_nmp*PwM2; % open loop function PI
L_NMP_Hinf = KHinf7*PwM2; % open loop function H-infinity
CL_NMP_PI = feedback(L_NMP_PI, eye(2)); % closed loop function PI
CL_NMP_Hinf= feedback(L_NMP_Hinf, eye(2)); % closed loop function H-infinity
[stabmarg1,wcu1] = robstab(CL_NMP_PI)
stabmarg1 =
LowerBound: 1.7651
UpperBound: 1.7687
CriticalFrequency: 0.0052
wcu1 =
Gs_nmp_unc_sampled_InputMultDelta: [2×2 ss]
[stabmarg2,wcu2] = robstab(CL_NMP_Hinf)
stabmarg2 =
LowerBound: 3.6340
UpperBound: 3.6413
CriticalFrequency: 8.9079e-04
wcu2 =
Gs_nmp_unc_sampled_InputMultDelta: [2×2 ss]
Being the upper and lower bounds striclty >1 in both case, we can state that both controllers are able to robust stabilize the uncertain system with MP property.
Let us check if they can also guarantee robust performance, i.d. if 
where M is the
component of the augmented plant
the only possible source of instability of the overall control structured system
and
is the structured singular value.
[perfmarg_MP_Hinf,wcu,report,info1] = robustperf(CL_NMP_Hinf)
perfmarg_MP_Hinf =
LowerBound: 0.6966
UpperBound: 0.6966
CriticalFrequency: 0.0647
wcu =
Gs_nmp_unc_sampled_InputMultDelta: [2×2 ss]
report =
'Uncertain system does not achieve performance robustness to modeled uncertainty. '
' -- The tradeoff of model uncertainty and system gain is balanced at a level of 69.7% of the modeled uncertainty. '
' -- A model uncertainty of 69.7% can lead to input/output gain of 1.44 at 0.0647 rad/seconds. '
' -- Sensitivity with respect to each uncertain element is: '
' 23% for Gs_nmp_unc_sampled_InputMultDelta. Increasing Gs_nmp_unc_sampled_InputMultDelta by 25% decreases the margin by 5.75%.'
info1 =
Sensitivity: [1×1 struct]
Frequency: [100×1 double]
BadUncertainValues: [100×1 struct]
MussvBnds: [1×2 frd]
MussvInfo: [1×1 struct]
[perfmarg_MP_PI,wcu,report,info2] = robustperf(CL_NMP_PI)
perfmarg_MP_PI =
LowerBound: 0.3708
UpperBound: 0.3709
CriticalFrequency: 0.0066
wcu =
Gs_nmp_unc_sampled_InputMultDelta: [2×2 ss]
report =
'Uncertain system does not achieve performance robustness to modeled uncertainty. '
' -- The tradeoff of model uncertainty and system gain is balanced at a level of 37.1% of the modeled uncertainty. '
' -- A model uncertainty of 37.1% can lead to input/output gain of 2.7 at 0.00658 rad/seconds. '
' -- Sensitivity with respect to each uncertain element is: '
' 20% for Gs_nmp_unc_sampled_InputMultDelta. Increasing Gs_nmp_unc_sampled_InputMultDelta by 25% decreases the margin by 5%.'
info2 =
Sensitivity: [1×1 struct]
Frequency: [87×1 double]
BadUncertainValues: [87×1 struct]
MussvBnds: [1×2 frd]
MussvInfo: [1×1 struct]
As in the MP case, while we had robust stability for both control systems, we don't have robust performance capabilities, being the lower and the upper bounds for robust performance strictly <1 in both cases.
Let us visualize how the structured singulare value varies with variation of frequencies.
% extracting the values for the H inf controller
a = info1.MussvBnds(1,1);
b = info1.MussvBnds(1,2);
max_nmp_Hinf = zeros(length(info1.Frequency),1);
for i=1:size(info1.Frequency)
max_nmp_Hinf(i)= a.Response(:,:,i);
min_nmp_Hinf = zeros(length(info1.Frequency),1);
for i=1:length(info1.Frequency)
min_nmp_Hinf(i)= b.Response(:,:,i);
% extracting the values for the PI controller
c = info2.MussvBnds(1,1);
d = info2.MussvBnds(1,2);
max_nmp_PI = zeros(length(info2.Frequency),1);
for i=1:length(info2.Frequency)
max_nmp_PI(i)= c.Response(:,:,i);
min_nmp_PI = zeros(length(info2.Frequency),1);
for i=1:length(info2.Frequency)
min_nmp_PI(i)= d.Response(:,:,i);
semilogx(info1.Frequency, max_nmp_Hinf,'b--',info1.Frequency, min_nmp_Hinf,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
title('µ bounds H-inf CL');
semilogx(info2.Frequency, max_nmp_PI,'b--',info2.Frequency, min_nmp_PI,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
Here we have another important limitation of the NMP case: the robust performance condition is much more violated by both controllers than the MP case. Anyaway the
again looks like to behave better than the PI. Let us try to design a
controller able to guarantee robust performance trhough μ-synthesis control. We have to build the augumented plant and run the D-K algorithm: the requirements on the weights have to be coherent with the limitations imposed by the NMP property of the system.
omb8 = 0.0064; % maximal bandiwith allowed for NMP zero
wp66 = ((s/M8+ omb8)/ (s+omb8*A8))* [0,1;1,0]; % off diagonal weight
% inizializing the inputs and the outputs of each block
PwM2.u = 'u'; % control input
PwM2.y = 'y'; % not perturbed output
wp66.y = 'z'; % performance variable
Sum4 = sumblk('v = r - y', 2); % measured variables
Pe_NMP = connect(PwM2,wp66,Sum4, {'r','u'},{'z','v' }); % augmented plant
ncont = 2; % two control inputs
nmeas = 2; % two measurements
[Krob2,rpMU2] = musyn(Pe_NMP ,nmeas,ncont)
D-K ITERATION SUMMARY:
-----------------------------------------------------------------
Robust performance Fit order
-----------------------------------------------------------------
Iter K Step Peak MU D Fit D
1 1.522 1.391 1.405 4
2 1.265 1.222 1.234 12
3 1.206 1.19 1.202 12
4 1.199 1.186 1.198 12
5 1.2 1.188 1.2 12
Best achieved robust performance: 1.19
Krob2 =
A =
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26
x1 -0.01107 0.23 0.1125 -0.006708 -0.08067 -0.04269 0.4043 -0.3527 -0.8182 0.3497 0.8875 -0.6259 -0.01516 -0.01827 0.01828 0.01098 0.01482 -0.0108 -0.4932 0.7911 0.004226 0.05411 0.02667 -0.001738 -0.02093 -0.01107
x2 -1.89e-16 -0.03478 1 1.753e-16 2.156e-15 1.137e-15 -1.462e-14 6.226e-15 1.774e-14 -6.438e-15 -2.581e-14 1.466e-14 2.434e-16 2.944e-16 -2.94e-16 -2.4e-16 -3.243e-16 2.339e-16 1.273e-14 -1.464e-14 -4.721e-17 -5.763e-16 -3.043e-16 4.399e-17 5.403e-16 2.849e-16
x3 0.07354 0.6774 -1.074 -0.09991 -1.202 -0.6359 6.022 -5.254 -12.19 5.209 13.22 -9.323 -0.2258 -0.2721 0.2723 0.1635 0.2208 -0.1608 -7.346 11.78 0.06295 0.8059 0.3972 -0.02588 -0.3117 -0.1649
x4 -0.006885 -0.07939 -0.04219 -0.01148 0.2242 0.1124 -0.6562 0.3644 0.741 -0.3956 -1.266 0.5324 0.01256 0.01527 -0.01519 -0.01657 -0.02239 0.01598 0.7574 -0.6578 -0.001834 -0.02137 -0.01134 0.004165 0.05318 0.02682
x5 1.029e-16 1.255e-15 6.625e-16 -9.538e-17 -0.03478 1 7.943e-15 -3.388e-15 -9.661e-15 3.504e-15 1.403e-14 -7.987e-15 -1.325e-16 -1.602e-16 1.6e-16 1.306e-16 1.764e-16 -1.273e-16 -6.925e-15 7.969e-15 2.57e-17 3.137e-16 1.657e-16 -2.393e-17 -2.939e-16 -1.55e-16
x6 -0.1026 -1.182 -0.6284 0.06744 0.591 -1.075 -9.773 5.428 11.04 -5.893 -18.86 7.93 0.1871 0.2274 -0.2263 -0.2469 -0.3335 0.2381 11.28 -9.799 -0.02732 -0.3182 -0.1689 0.06203 0.7922 0.3996
x7 0.04902 0.4652 0.5305 -0.06672 -0.8049 -0.4257 4.136 -3.353 -8.043 3.34 8.858 -6.242 -0.1401 -0.1693 0.1692 0.1068 0.1442 -0.1049 -4.895 7.583 0.01464 0.1512 0.1398 -0.01721 -0.2078 -0.1099
x8 0.0756 0.6569 0.9166 -0.1127 -1.354 -0.7167 6.576 -5.996 -13.99 5.915 14.66 -10.74 -0.2512 -0.3037 0.3036 0.1851 0.2499 -0.182 -8.282 13.43 0.02367 0.234 0.2435 -0.02921 -0.3516 -0.186
x9 0.02788 0.3401 0.1796 -0.02583 -0.3176 -0.1674 2.147 -0.9025 -2.62 0.9493 3.794 -2.167 -0.0359 -0.04342 0.04336 0.03538 0.0478 -0.03448 -1.876 2.16 0.006966 0.08503 0.0449 -0.006481 -0.07961 -0.04198
x10 -0.114 -1.312 -0.6976 0.06679 0.5354 0.9045 -11.1 6.103 12.26 -6.671 -21.38 8.75 0.2099 0.2552 -0.2539 -0.2728 -0.3689 0.2638 12.72 -10.95 -0.03043 -0.354 -0.188 0.02215 0.2124 0.2431
x11 -0.01345 -0.164 -0.0866 0.01249 0.1536 0.08097 -1.044 0.4432 1.261 -0.4427 -1.842 1.042 0.01733 0.02096 -0.02093 -0.0171 -0.0231 0.01666 0.9069 -1.042 -0.00336 -0.04101 -0.02166 0.003134 0.03849 0.0203
x12 -0.07143 -0.832 -0.4416 0.04678 0.4323 0.539 -6.667 3.527 7.481 -3.806 -12.65 5.481 0.1235 0.15 -0.1494 -0.1547 -0.2092 0.1498 7.331 -6.589 -0.01882 -0.221 -0.1172 0.01443 0.1475 0.1433
x13 -0.06701 -0.838 -0.4262 0.04465 0.5466 0.2884 -3.979 1.747 4.154 -1.777 -7.049 3.212 -0.05918 0.2008 -0.1379 -0.06365 -0.08597 0.06217 3.263 -3.886 -0.1765 -2.485 -0.8002 0.01128 0.138 0.07282
x14 0.05024 0.6232 0.3209 -0.0378 -0.4635 -0.2445 3.269 -1.424 -3.653 1.458 5.789 -2.911 0.01952 -0.1494 0.09497 0.05302 0.07162 -0.05175 -2.755 3.24 0.09243 1.294 0.4271 -0.009522 -0.1167 -0.06156
x15 -0.01066 -0.1326 -0.06798 0.007676 0.09406 0.04962 -0.6681 0.2932 0.7355 -0.2995 -1.185 0.581 -0.007333 0.0401 -0.0437 -0.01083 -0.01463 0.01057 0.56 -0.6627 -0.02263 -0.3176 -0.1037 0.001936 0.0237 0.01251
x16 0.04815 0.5784 0.3059 -0.06938 -0.8721 -0.4477 3.399 -1.867 -4.919 1.975 6.484 -3.958 -0.0693 -0.08397 0.08374 -0.07588 0.2361 -0.1409 -3.742 4.022 0.01226 0.1476 0.07805 -0.2025 -2.856 -0.9161
x17 -0.03762 -0.4544 -0.2402 0.0473 0.5911 0.3056 -2.758 1.38 3.727 -1.45 -5.107 3.025 0.05211 0.06311 -0.06296 0.03215 -0.1643 0.08726 2.787 -3.06 -0.009513 -0.1151 -0.06085 0.1044 1.467 0.4785
x18 0.00845 0.1019 0.05389 -0.01095 -0.137 -0.07071 0.6177 -0.3137 -0.8414 0.33 1.149 -0.6813 -0.0118 -0.01429 0.01426 -0.009242 0.04517 -0.04356 -0.6329 0.6909 0.00214 0.02587 0.01367 -0.02589 -0.3641 -0.1182
x19 -2.435e-13 -2.97e-12 -1.568e-12 2.258e-13 2.776e-12 1.464e-12 -1.881e-11 8.021e-12 2.287e-11 -8.294e-12 -3.321e-11 1.89e-11 3.136e-13 3.793e-13 -3.788e-13 -3.092e-13 -4.177e-13 3.014e-13 -0.00064 -1.886e-11 -6.084e-14 -7.427e-13 -3.922e-13 5.665e-14 6.959e-13 3.67e-13
x20 6.803e-15 8.298e-14 4.381e-14 -6.308e-15 -7.756e-14 -4.089e-14 5.255e-13 -2.241e-13 -6.388e-13 2.317e-13 9.278e-13 -5.281e-13 -8.761e-15 -1.06e-14 1.058e-14 8.638e-15 1.167e-14 -8.419e-15 -4.58e-13 -0.00064 1.7e-15 2.075e-14 1.096e-14 -1.583e-15 -1.944e-14 -1.025e-14
x21 0.3046 3.714 1.962 -0.2835 -3.485 -1.838 23.52 -10.07 -28.81 10.41 41.57 -23.85 -0.3932 -0.4757 0.4749 0.3881 0.5243 -0.3782 -20.58 23.7 0.04716 0.929 0.4908 -0.07112 -0.8736 -0.4607
x22 0.02298 0.2808 0.1482 -0.02009 -0.2466 -0.13 1.459 -0.741 -2.272 0.7633 2.698 -1.922 -0.02913 -0.03523 0.03518 0.02792 0.03772 -0.02723 -1.46 1.789 0.005727 0.03525 1.037 -0.005053 -0.06197 -0.03269
x23 -0.07915 -0.9912 -0.5037 0.05042 0.6162 0.3252 -4.174 2.027 5.038 -2.057 -7.588 3.965 0.08472 0.1023 -0.1023 -0.0727 -0.09819 0.07106 3.688 -4.574 -0.2122 -2.989 -2.872 0.01276 0.1559 0.0823
x24 -0.2405 -2.934 -1.549 0.2218 2.725 1.438 -18.72 7.91 22.52 -8.177 -32.97 18.61 0.3094 0.3743 -0.3737 -0.3046 -0.4115 0.2969 16.19 -18.61 -0.06009 -0.7336 -0.3873 0.0267 0.6835 0.3607
x25 -0.02375 -0.2893 -0.1528 0.02282 0.2808 0.148 -2.04 0.7928 2.156 -0.8216 -3.523 1.754 0.03088 0.03737 -0.0373 -0.03098 -0.04185 0.03018 1.656 -1.834 -0.005943 -0.07247 -0.03827 0.005717 0.0355 1.037
x26 0.0473 0.5672 0.3 -0.07072 -0.8897 -0.4565 3.62 -1.865 -4.746 1.978 6.777 -3.77 -0.06893 -0.08354 0.0833 0.08063 0.1089 -0.07818 -3.787 3.959 0.01207 0.1451 0.07674 -0.2102 -2.965 -2.861
B =
u1 u2
x1 -1.656e-14 2.029e-14
x2 -1.842e-14 2.233e-14
x3 -2.539e-15 3.117e-15
x4 -6.646e-15 6.774e-15
x5 9.966e-15 -1.221e-14
x6 2.487e-16 -3.136e-16
x7 0.5508 -1.76
x8 -1.054 -0.8439
x9 2.68 -3.325
x10 -0.3871 -0.8198
x11 -1.325 1.578
x12 -1.613 1.283
x13 -6.085 0.9564
x14 4.616 -2.502
x15 -0.9605 0.4071
x16 0.8523 -5.454
x17 -1.788 4.296
x18 0.3617 -0.9548
x19 0.0625 2.888e-11
x20 6.597e-13 0.0625
x21 29.19 -36.99
x22 1.016 -3.683
x23 -5.263 2.481
x24 -24.03 28.2
x25 -3.093 2.192
x26 1.767 -4.435
C =
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26
y1 0.3471 3.197 3.947 -0.4715 -5.671 -3.001 28.42 -24.79 -57.51 24.58 62.38 -44 -1.065 -1.284 1.285 0.7716 1.042 -0.7589 -34.67 55.61 0.2971 3.803 1.875 -0.1221 -1.471 -0.7784
y2 -0.484 -5.58 -2.966 0.3183 2.789 3.947 -46.12 25.62 52.09 -27.81 -89 37.43 0.8829 1.073 -1.068 -1.165 -1.574 1.123 53.24 -46.24 -0.1289 -1.502 -0.7973 0.2927 3.738 1.886
D =
u1 u2
y1 0 0
y2 0 0
Continuous-time state-space model.
Model Properties
rpMU2 = 1.1862
After only 5 iterations, the best attempt of musyn() using D-K algorithm ends with a
order controller able to guarantee a maximum structured singular value of 1.18, so not robust performant. Also the attempt of relaxing the requirements in the weight selections parameters is not sufficient to achieve better performances: we accepted also 1% of steady state error and an maximal overshoot of 2. We expect lower and upper bounds for robust performance indices significantly lower than the MP case.
L_NMP_rob = Krob2*PwM2; % open loop function robust performance controller (MP)
CL_NMP_rob = feedback(L_NMP_rob, eye(2)); % closed loop robust performance controller (MP)
[perfmarg_NMP_rob,wcu,report,info4] = robustperf(CL_NMP_rob)
perfmarg_NMP_rob =
LowerBound: 0.7627
UpperBound: 0.7627
CriticalFrequency: 0.0875
wcu =
Gs_nmp_unc_sampled_InputMultDelta: [2×2 ss]
report =
'Uncertain system does not achieve performance robustness to modeled uncertainty. '
' -- The tradeoff of model uncertainty and system gain is balanced at a level of 76.3% of the modeled uncertainty. '
' -- A model uncertainty of 76.3% can lead to input/output gain of 1.31 at 0.0875 rad/seconds. '
' -- Sensitivity with respect to each uncertain element is: '
' 24% for Gs_nmp_unc_sampled_InputMultDelta. Increasing Gs_nmp_unc_sampled_InputMultDelta by 25% decreases the margin by 6%.'
info4 =
Sensitivity: [1×1 struct]
Frequency: [89×1 double]
BadUncertainValues: [89×1 struct]
MussvBnds: [1×2 frd]
MussvInfo: [1×1 struct]
Let us visualize the behavior of the structured singular value in a suitable range of frequencies.
% extracting the structured singular values for the µ-syn robust controller
p = info4.MussvBnds(1,1);
o = info4.MussvBnds(1,2);
max_nmp_rob = zeros(length(info4.Frequency),1);
for i=1:length(info4.Frequency)
max_nmp_rob(i)= p.Response(:,:,i);
min_nmp_rob = zeros(length(info4.Frequency),1);
for i=1:length(info4.Frequency)
min_nmp_rob(i)= o.Response(:,:,i);
semilogx(info4.Frequency, max_nmp_rob,'b--',info4.Frequency, min_nmp_rob,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
title('µ bounds robust performant CL');
We can see that µ is not always <1: in a significantly large neighborood of the critical frequency
the structured singular value is about 1,27. This result is reasonably the best achieved in terms of robust performance, better than the
nominal controller which achieves a robust performance index around 69.66%, and than the PI that achieves about 37,08%, but it is consistently worse than the MP case. This is another symptom of how much difficult is to control a NMP MIMO system than one will all zeros 
Let us compute the step response and try to visualize some possible behaviors (10) of closed loop system under the µ-syn controller.
L_NMP_rob = Krob2* PwM2; % robust open loop
CL_NMP_rob = feedback(L_NMP_rob, eye(2)); % robust closed loop
r_unit4 = [ones(length(t1),1),ones(length(t1),1)];
samples_NMP = usample(CL_NMP_rob ,10);
lsim(samples_NMP , 'b',r_unit4, t1);
title("System's steady state with µ-synthesis controller (NMP)");
In the end, we can note two intresting aspects:
- as expected, and resonable, from the plots it is much more evident the uncertainty nature of the closed loop system, then the MP case, because the system samples tends to assume slightly different behaviors, thanks to the lower robust performance margins;
- althought this, the system will never diverge under such proposed feedbacks since all of them are proved to be robust stable;
- the closed-loop outputs under the controller designed with the µ-synthesis exhibits a transient faster than the case of PI and the
nominal controller. This is quite surprising since the required bandiwith in the design of the performance weight is the same imposed by the NMP zero's limitation.
This concludes our robustness analysis.
5. Decoupling approaches.
Let us try to design two possible decoupler in order to separate the input-output channels
minimizing or hopefully nihilating the interactions. We want to find the matrix
such that the MIMO system behaves like two decoupled SISO systems, with the possibility of controlling them with two additional SISO controllers (Simplified Decoupling approach), or including the control in the design of
(Inverted Decoupling approach). In the first case the overall control system will look like:
while in the second one
will substantially coincide with 
Simplified Decoupling
The first attempt is based on what it is called "Simplified feedforward decoupling". Matrix
is composed in this way: and the overall "apparent process" will be:
. Note that the zeros of
are the same as
. Furthermore, since both
and
are non-zero, all transmition zeros not coinciding with poles (i.e. solutions of
, will appear in the resulting open loop system as zeros of the elements on the diagonal. This is an important limitation, since we have a case in which the system's transfer function exhibitis a Non Minimum Phase transmition zero.
So we expect two positive zeros, one on each element of the diagonal, that are more difficult to manage since NMP zeros in the SISO case don't have directions and their effects are present in all situations, while in the MIMO case they act only along the zeros directions.
d12mp = - (Gs_mp(1,2)/(Gs_mp(1,1))); %element out of the diagonal of D(s)
d21mp= - (Gs_mp(2,1)/(Gs_mp(2,2))); %element out of the diagonal of D(s)
Ds_mp = [1, d12mp; d21mp, 1];
L_mp_app= simplify(zpk(Gs_mp*Ds_mp)) % open loop
L_mp_app =
From input 1 to output...
0.041625 (s+0.05802) (s+0.01718) (s+0.01595) (s+0.01107)
1: --------------------------------------------------------
(s+0.04186) (s+0.03334) (s+0.01595)^2 (s+0.01107)
6.1876e-20 s (s+0.07544) (s-0.01294)
2: ------------------------------------
(s+0.03334)^2 (s+0.01107)^3
From input 2 to output...
-1.585e-19 s (s^2 + 0.001953)
1: -----------------------------
(s+0.04186)^2 (s+0.01595)^3
0.031406 (s+0.05802) (s+0.01718) (s+0.01595) (s+0.01107)
2: --------------------------------------------------------
(s+0.04186) (s+0.03334) (s+0.01595) (s+0.01107)^2
Continuous-time zero/pole/gain model.
Model Properties
Since the elements off-diagonal exhibits coefficient of the order of
, we can consider them identically equal to zero (they are actual different just for numerical issues for MATLAB). So we can consider our decoupling problem solved. In order to avoid numerical problems in the next few computations, let's substitute zero off-diagonal instead of the remaining very small values, whose presence is a consequence of numerical MATLAB approximations.
L_mp= [L_mp_app(1,1), 0; 0, L_mp_app(2,2)]; % open looop actually decoupled
zero(L_mp(1,1))
-0.0580
-0.0172
-0.0159
-0.0111
zero(L_mp(2,2))
-0.0580
-0.0172
-0.0159
-0.0111
pole(L_mp(1,1))
-0.0419
-0.0333
-0.0159
-0.0159
-0.0111
pole(L_mp(2,2))
-0.0419 + 0.0000i
-0.0333 + 0.0000i
-0.0159 + 0.0000i
-0.0111 + 0.0000i
-0.0111 - 0.0000i
We can see from the poles and the zeros computation that there are several cancellations in the remaining two SISO elements on the diagonal, but this fact does not constitute a problem since they all are with negative real part. The design of a stabilizing controller can start from here.
Let us explore the Non Minimum Phase Case.
d12nmp = - (Gs_nmp(1,2)/(Gs_nmp(1,1)));
d21nmp= - (Gs_nmp(2,1)/(Gs_nmp(2,2)));
Ds_nmp = zpk([1, d12nmp; d21nmp, 1])
Ds_nmp =
From input 1 to output...
1: 1
-0.028515 (s+0.01094)
2: -----------------------
(s+0.01782) (s+0.01094)
From input 2 to output...
-0.041223 (s+0.01582)
1: -----------------------
(s+0.02563) (s+0.01582)
2: 1
Continuous-time zero/pole/gain model.
Model Properties
L_nmp_app = simplify(zpk(Gs_nmp*Ds_nmp))
L_nmp_app =
From input 1 to output...
0.024111 (s+0.05623) (s+0.01582) (s+0.01094) (s-0.01278)
1: --------------------------------------------------------
(s+0.01582)^2 (s+0.01782) (s+0.02563) (s+0.01094)
1.4677e-19 (s-0.000503) (s+0.01094) (s+0.02335)
2: -----------------------------------------------
(s+0.01782)^2 (s+0.01094)^3
From input 2 to output...
-3.5894e-19 (s+0.01582) (s^2 + 0.0215s + 0.0002711)
1: ---------------------------------------------------
(s+0.01582)^3 (s+0.02563)^2
0.017478 (s+0.05623) (s+0.01582) (s+0.01094) (s-0.01278)
2: --------------------------------------------------------
(s+0.01782) (s+0.01582) (s+0.02563) (s+0.01094)^2
Continuous-time zero/pole/gain model.
Model Properties
Again we have the same numerical problems as before.
L_nmp= [L_nmp_app(1,1), 0; 0, L_nmp_app(2,2)]; % open looop actually decoupled
zero(L_nmp(1,1))
-0.0562
0.0128
-0.0158
-0.0109
zero(L_nmp(2,2))
-0.0562
0.0128
-0.0109
-0.0158
pole(L_nmp)
-0.0158
-0.0256
-0.0158
-0.0178
-0.0109
-0.0178
-0.0109
-0.0109
-0.0256
-0.0158
We can see the presene of NMP zeros in the two transfer functions, as well as the presence of several cancellations.
Inverted Decoupling
The goal now is to design a controller able to solve the decoupling problem and at the same time to impose a desired behavior.
So let us choose a possible desired open loop to impose through a positive inner feedback (2x2) decoupler of the form
following the control scheme:
The choices
,
,
,
, guarantee that
, with
the desired behavior of the open loop, decoupled and stabilizated. The resulting control block scheme is:
We choose the desired open loop such that it is diagonal, type 1, with all poles
i.e. close loop stable and corresponding to two equivalent and indipendent SISO systems. Its transfer function is:

Ldes = [100/(s*(s+20)), 0 ;0 ,100/(s*(s+20))]; % inizializing the desired open loop transfer function
Dd11 = Ldes(1,1)/ Gs_nmp(1,1);
Dd22 = Ldes(2,2)/ Gs_nmp(2,2);
Do12 = - (Gs_nmp(1,2)/(Ldes(1,1)));
Do21 = - (Gs_nmp(2,1)/( Ldes(2,2)));
Dd = [Dd11, 0;0, Dd22]; % building Dd(s)
Do = [0, Do12; Do21, 0]; % building Do(s)
D_ID = Dd*(inv(eye(2)-(Do*Dd))); % building D(s)
L_ID_app = Gs_nmp * D_ID % computing G(s) D(s)
L_ID_app =
From input 1 to output...
3.794e07 s^17 + 3.802e09 s^16 + 1.525e11 s^15 + 3.065e12 s^14 + 3.096e13 s^13 + 1.275e14 s^12 + 2.447e13 s^11 + 1.871e12 s^10 + 6.972e10 s^9 + 1.16e09 s^8 - 1.528e06 s^7 - 3.709e05 s^6 - 5055 s^5 - 4.737 s^4 + 0.4714 s^3 + 0.004267 s^2 + 1.206e-05 s
1: ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
3.794e05 s^19 + 4.561e07 s^18 + 2.286e09 s^17 + 6.116e10 s^16 + 9.227e11 s^15 + 7.466e12 s^14 + 2.574e13 s^13 + 4.913e12 s^12 + 3.75e11 s^11 + 1.396e10 s^10 + 2.32e08 s^9 - 3.094e05 s^8 - 7.423e04 s^7 - 1011 s^6 - 0.9426 s^5 + 0.09432 s^4 + 0.0008536 s^3 + 2.411e-06 s^2
1.546e-06 s^16 + 0.0001238 s^15 + 0.003719 s^14 + 0.04965 s^13 + 0.2498 s^12 + 0.008789 s^11 - 0.003662 s^10 - 0.0004425 s^9 - 2.05e-05 s^8 - 4.694e-07 s^7 - 4.54e-09 s^6 + 3.638e-11 s^5 + 1.464e-12 s^4 + 1.643e-14 s^3 + 8.674e-17 s^2 + 1.813e-19 s
2: ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1.141e06 s^19 + 1.371e08 s^18 + 6.871e09 s^17 + 1.838e11 s^16 + 2.771e12 s^15 + 2.24e13 s^14 + 7.701e13 s^13 + 1.342e13 s^12 + 9.03e11 s^11 + 2.798e10 s^10 + 3.234e08 s^9 - 3.443e06 s^8 - 1.406e05 s^7 - 1258 s^6 + 4.036 s^5 + 0.1511 s^4 + 0.001031 s^3 + 2.411e-06 s^2
From input 2 to output...
-2.375e-07 s^16 - 1.916e-05 s^15 - 0.0005817 s^14 - 0.007965 s^13 - 0.04272 s^12 - 0.02441 s^11 - 0.003906 s^10 - 0.0003052 s^9 - 1.431e-05 s^8 - 3.278e-07 s^7 - 3.725e-09 s^6 - 2.728e-12 s^5 + 4.974e-13 s^4 + 6.439e-15 s^3 + 3.469e-17 s^2 + 7.962e-20 s
1: ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
3.794e05 s^19 + 4.561e07 s^18 + 2.286e09 s^17 + 6.116e10 s^16 + 9.227e11 s^15 + 7.466e12 s^14 + 2.574e13 s^13 + 4.913e12 s^12 + 3.75e11 s^11 + 1.396e10 s^10 + 2.32e08 s^9 - 3.094e05 s^8 - 7.423e04 s^7 - 1011 s^6 - 0.9426 s^5 + 0.09432 s^4 + 0.0008536 s^3 + 2.411e-06 s^2
1.141e08 s^17 + 1.143e10 s^16 + 4.585e11 s^15 + 9.211e12 s^14 + 9.293e13 s^13 + 3.817e14 s^12 + 6.69e13 s^11 + 4.508e12 s^10 + 1.398e11 s^9 + 1.618e09 s^8 - 1.718e07 s^7 - 7.027e05 s^6 - 6291 s^5 + 20.14 s^4 + 0.7551 s^3 + 0.005153 s^2 + 1.206e-05 s
2: ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1.141e06 s^19 + 1.371e08 s^18 + 6.871e09 s^17 + 1.838e11 s^16 + 2.771e12 s^15 + 2.24e13 s^14 + 7.701e13 s^13 + 1.342e13 s^12 + 9.03e11 s^11 + 2.798e10 s^10 + 3.234e08 s^9 - 3.443e06 s^8 - 1.406e05 s^7 - 1258 s^6 + 4.036 s^5 + 0.1511 s^4 + 0.001031 s^3 + 2.411e-06 s^2
Continuous-time transfer function.
Model Properties
Since we have the usual numerical problems, we substitute zero at the off-diagonal positions.
Let us now compare the resulting open loop with the desired one from a frequency domain point of view.
L_ID= [L_ID_app(1,1), 0;0, L_ID_app(2,2)]; % resultingo L(s) with 0 off-diagonal
Tit = title('Magnitude of Desired Open Loop');
set(Tit,'fontName','Courier','FontSize',10);
Tit = title('Magnitude of Obtained Open Loop');
set(Tit,'fontName','Courier','FontSize',10);
phdes = bodeplot(Ldes); grid on
setoptions(phdes,'MagVisible', 'off');
Tit = title('Phase of Desired Open Loop');
set(Tit,'fontName','Courier','FontSize',10);
phobt = bodeplot(L_ID); grid on
setoptions(phobt,'MagVisible', 'off');
Tit = title('Phase of Obtained Open Loop');
set(Tit,'fontName','Courier','FontSize',10);
First thing we note is that while the two SISO model on the diagonal of the obtained open loop exhibits the same phase, they differ a bit from the point of view of the magnitude, signal that we are made little numerical errors in computations.
The form of the phases and of the magnitudes are more or less the same as the ones of the desired open loops, but they result a bit anticipated in frequency, exept for the magnitude of the second element on the diagonal, that it exactly the same as the desired one.
Let us evaluate the performances in terms of stability and reference reproduction's capabilities.
T_ID = L_ID(1,1)/(1+L_ID(1,1)); % complementary sensitivity
r_unit3 = [ones(length(t),1)]; % dummy step reference to check if the controller works
ycl3 = lsim(T_ID, r_unit3, t);
plot(t, ycl3, 'LineWidth', 2), grid on;
title("Steady state for each decoupled I/O channel");
legend('y1(t)= y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
So as expected the inverted decoupler is able to make the closed loop system decoupled and exhibiting a behavior close to the desired, although the overall open loop is not excactly the same as the desired one.
We expect an high control effort peak due to the fact that we are decoupling and stabilizing at the same time.
Cso = D_ID/(eye(2)+ L_ID_app); % control sensitivity function
m0 = lsim(Cso, r_unit2, t);
plot(t,m0,'LineWidth',2); grid
title('Control input for r_{i}(t) = step with Inverted Decoupler','FontName','courier','FontSize',12)
xlabel('time','FontName','courier','FontSize',14);
ylabel('u_{i}(t)','FontName','courier','FontSize',14);
legend('u1(t)',' u2(t)','FontName','courier','FontSize',12,'Location','NorthEast')